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ABSTRACT 


This  thesis  dev  an  orbital  prediction  model  based  on  fundamental  principles  of 
orbital  dynamics  and  drag.  A  FORTRAN  based  orbital  prediction  scheme  was  designed  to 
provide  accurate  ephemerides  for  a  particular  DoD  satellite  program.  The  satellite 
program  under  study  has  satellites  at  650  and  800  kilometers  with  high  inclinations.  In 
order  to  obtain  the  highest  accuracy  possible,  a  comparison  of  atmospheric  models  had  to 
be  conducted  in  order  to  determine  which  model  was  more  accurate.  Mathematical 
formulation  for  three  widely  used  earth  atmospheric  models  are  presented;  the  JACCHIA 
60,  JACCHIA  71,  and  MSIS  86  atmosphe  ■'  v  s.  The  MSIS  86  atmospheric  model 
was  not  evaluated  due  to  computer  problems.  C  omptrison  of  the  two  JACCHIA  models 
proved  that  the  JACCHIA  71  model  provided  much  more  accurate  ephemerides.  It  is 
believed  that  this  is  due  not  only  to  the  incorporation  of  variations  in  densi.y  caused  by 
solar  flux,  but  also  geomagnetic  activity  and  a  better  modeling  of  the  polar  regions. 

Further  work  on  this  project  would  include  incorporation  of  the  MSIS  86  model  for 
evaluation,  incorporation  of  the  full  WGS-84  geopotential  model,  and  using  more  accurate 
observed  vectors  in  order  to  obtain  a  better  comparison.  Incorporating  a  subroutine 
which  will  vary  the  B-factor  as  a  function  of  latitude  will  greatly  increase  accuracy.  This  is 
a  major  deviation  from  current  operational  practice,  in  that  the  B-factor  is  often  used  as  an 
error  catch-all  and  does  not  truly  represent  its  dynamical  purpose. 
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I.  INTRODUCTION 


Orbital  prediction  has  become  an  essential  science  needed  for  several  of  the  DoD 
satellite  programs.  Exact  satellite  ephemerides  provide  for  a  more  accurate  means  of 
mission  analysis.  Put  more  directly,  the  more  accurate  the  satellite  ephemerides  can  be 
calculated  or  predicted,  the  more  accurate  the  mission  analysis  becomes.  The  main 
objective  being  to  pinpoint  the  satellites  current  or  future  position.  Since  the  launch  of 
Sputnik  in  the  late  fifties,  a  great  deal  of  effort  has  been  placed  in  trying  to  model  the 
space  environment.  In  particular,  the  upper  atmosphere,  or  neutral  atmosphere,  has  been 
of  great  interest  in  the  study  of  artificial  satellite  orbit  theory.  Since  1957,  several 
attempts  have  been  made  in  modeling  the  thermosphere  in  order  to  aid  in  satellite  mission 
analysis.  Several  atmospheric  models,  or  satellite  drag  models,  are  currently  used  for 
practical  applications  such  as  lifetime  estimates,  reentry  prediction,  orbit  determination  and 
tracking,  attitude  dynamics,  and  most  recently,  mass  analysis  of  a  particular  satellite. 
Atmospheric  drag  affects  all  satellites,  in  all  altitude  regimes,  from  low  earth  orbits  to  well 
beyond  geosynchronous  altitudes.  For  many  satellites,  the  modeling  of  atmospheric  drag 
is  the  largest  error  source  in  describing  the  forces  acting  on  the  satellite. 

Satellite  drag  models  can  be  divided  into  two  categories,  the  empirical  models  and 
the  general  circulation  models.  This  thesis  will  compare  three  of  the  more  widely  used 
empirical  models.  The  models  used  for  the  comparison  are  the  Jacchia  60,  the  Jacchia  71, 
and  the  MSIS  86  earth  atmospheric  models.  Each  of  these  models  is  formulated  in  a 
different  manner  and  are  unique  for  the  altitude  bands  that  were  tested.  The  ultimate  goal 
being  to  find  out  which  model  is  more  accurate  in  terms  of  orbit  prediction  for  a  particular 
satellite  program.  Since  the  modeling  of  atmospheric  drag  is  the  largest  error  in  orbital 
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prediction,  finding  the  more  accurate  model  for  the  satellite  program  in  question  will 
greatly  improve  the  predicted  satellite  ephemirides  and  hence  mission  analysis 

The  atmospheric  models  being  used  in  this  comparison  ail  differ  with  respect  to  the 
input  that  each  model  requires.  Currently  the  Air  Force  Satellite  Tracking  Center  (STC) 
uses  the  Jacchia  60  atmospheric  model  which  uses  not  only  date  and  time  as  inputs  but 
also  the  measurement  of  the  solar  flux,  FI  0.7.  The  J71  model,  considered  to  be  an 
improvement  of  the  Jacchia  60  model,  adds  the  measurement  of  geomagnetic  activity,  or 
Ap,  to  its  inputs.  The  MSIS  86  atmospheric  model  is  formulated  in  a  different  manner 
than  the  Jacchia  atmospheric  models,  and  uses  both  FI 0.7,  and  Ap  as  its  inputs. 

The  atmospheric  models  calculate  the  density  and  constituency  of  the  atmosphere 
based  on  the  current  and  predicted  or  average  environmental  conditions.  The  accuracy  of 
these  models  has  been  calculated  to  be  80  to  85  percent  accurate.  This  percentage  drops 
off  substantially  as  the  altitude  increases.  Accuracy  also  seems  to  decrease  during  periods 
of  high  solar  and  geomagnetic  activity.  At  the  moment,  the  sun  is  on  the  downside  of  the 
1 1  year  solar  cycle.  This  is  an  advantage  for  the  comparison  of  the  atmospheric  models, 
because  there  were  fewer  and  weaker  solar  flares  and  geomagnetic  storms  to  perturb  the 
upper  atmosphere. 
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II.  BACKGROUND 


A.  ATMOSPHERE 

The  earth's  atmosphere  is  classically  divided  into  four  different  regions  based  on 
temperature  and  pressure  gradients.  These  four  regions  are  the  troposphere,  stratosphere, 
mesosphere  and  the  thermosphere.  The  corresponding  boundary  layers,  or  upper  limits  of 
each  of  the  regions  are  the  tropopause,  stratopause,  mesopause,  and  thermopause 
respectively.  Figure  1  (Akasofii,  1972,  p.  109),  illustrates  the  breakdown  of  the  earth's 
atmospheric  regions.  Beyond  the  thermopause  is  the  region  delineated  as  the  exosphere. 
The  exosphere  is  a  region  of  extremely  low  density  and  temperature,  and  is  the  transition 
region  into  space. 

Another  way  in  which  the  earth's  atmosphere  is  divided,  is  by  the  classification  of 
two  regions  known  as  the  homosphere  and  the  heterosphere.  The  transition  boundary 
between  the  regions  is  labeled  the  turbopause.  Once  again  the  region  outside  the 
heterosphere  is  labeled  the  exosphere.  Figure  2  (Akasofii,  1972,  p.  1 1 1)  illustrates  the 
various  atmospheric  regions  as  derived  by  the  two  defining  systems.  This  figure  provides 
a  breakdown  of  altitude  versus  temperature  for  the  various  altitude  regimes.  It  should  be 
noted  at  this  time,  that  the  atmospheric  region  of  interest  during  this  study  was  an  altitude 
band  between  600  to  800  km.  This  altitude  band  is  encompassed  within  either  the 
thermosphere  or  exosphere,  depending  on  which  classification  scheme  is  being  used.  For 
the  sake  of  continuity,  the  altitude  band  of  interest  will  be  considered  to  be  within  the 
thermosphere. 
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1.  Troposphere 

The  troposphere  is  that  portion  of  the  atmosphere  that  extends  from  the  surface 
to  roughly  10  to  15  km  above  the  surface.  It  is  in  equilibrium  with  the  sun-warmed 
surface  and  is  characterized  by  intense  convection  and  cloud  formations.  In  this  region, 
both  temperature  and  density  decrease  with  increasing  altitude,  with  an  occasional 
inversion  layer.  (U  S.  Air  Force,  1960,  p.  1-3) 

2.  Stratosphere 

The  stratosphere  is  located  above  the  tropopause  and  extends  up  to  50  km.  This 
region  of  the  atmosphere  is  extremely  important  in  that  it  contains  the  ozone  that  is 
responsible  for  the  absorption  of  the  extreme  ultraviolet  (EUV)  radiation  produced  by  the 
sun.  Due  to  the  absorption  of  this  EUV  radiation,  the  stratosphere  has  a  positive 
temperature  gradient.  The  density,  however,  still  decreases  with  altitude.  One  other 
consequence  of  the  absorption  of  the  EUV  radiation  by  the  ozone  layer  is  that  the  EUV 
radiation  cannot  be  measured  from  the  earth's  surface.  This  presents  a  problem  which  will 
be  addressed  later  in  this  paper. 

3.  Mesosphere 

The  Mesosphere  is  that  region  of  the  atmosphere  located  above  the  stratopause 
and  extends  up  to  80km.  Once  again,  both  the  temperature  and  density  are  decreasing 
with  altitude.  The  mesosphere  is  in  radiative  equilibrium  between  the  ultraviolet  ozone 
heating  by  the  upper  fringe  of  the  stratosphere,  and  the  infrared  ozone  and  carbon  dioxide 
cooling  by  radiation  to  space.  (U.S.  Air  Force,  1960,  p.  1-3) 

4.  Thermosphere 

The  thermosphere  extends  from  the  mesopause  to  higher  altitudes  with  no 
altitude  limit.  The  thermosphere  is  characterized  by  a  very  rapid  increase  in  temperature 
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4.  Thermosphere 

The  thermosphere  extends  from  the  mesopause  to  higher  altitudes  with  no 
altitude  limit.  The  thermosphere  is  characterized  by  a  very  rapid  increase  in  temperature 
with  altitude  due  to  the  absorption  of  the  sun's  EUV  radiation.  The  temperature  increase 
reaches  a  limiting  value  known  as  the  exospheric  temperature,  the  average  values  being 
between  -600  to  1200  K  over  a  solar  cycle.  (Larson,  1992,  p.208)  The  thermosphere  may 
also  be  heated  by  geomagnetic  activity,  which  transfers  energy  from  the  magnetosphere 
and  ionosphere  to  th*  ‘henmosphere.  The  heating  of  the  thermosphere  causes  an  increase 
in  the  atmospheric  density  due  to  the  expansion  of  the  atmosphere.  Figure  3  (Hess,  1965, 
p.  679),  illustrates  the  variation  of  temperature  versus  altitude  for  the  various  atmospheric 
regimes. 

5.  Homosphere 

The  homosphere  extends  from  the  surface  to  approximately  100km.  It  is 
characterized  by  its  uniform  composition  and  relatively  constant  molecular  weight.  The 
composition  of  this  region  can  be  broken  down  into  the  following:  78%  N,  21%  02,  1% 
At  and  trace  amounts  of  other  gases.  The  uniformity  of  the  region  is  created  due  to  the 
turbulent  mixing  of  the  gas  constituents.  (Adler,  1993,  p.  10)  The  composition,  hence  the 
uniformity,  of  the  homosphere  changes  at  -100km  altitude  due  mainly  to  the  dissociation 
of  the  oxygen  molecules.  Because  the  density  at  this  altitude  is  low,  recombination  of  the 
monatomic  oxygen  is  very  infrequent;  even  more  so  as  altitude  increases.  The  dissociation 
of  oxygen  causes  the  molecular  weight  to  decrease  substantially. 

6.  Heterosphere 

The  heterospheie  exists  from  ~  100km  outward,  with  no  altitude  limit.  The 
region  is  characterized  by  diffusive  equilibrium  and  significantly  varying  composition.  The 
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molecular  weight  of  the  atmosphere  decreases  rapidly,  from  ~  29  at  90km  to  -  16  at 
500km.  Above  the  region  of  oxygen  dissociation,  nitrogen  begins  to  dissociate.  Diffusive 
equilibrium  begins  to  take  place,  and  the  lighter  molecules  and  atoms  rise  to  the  top  of  the 
atmosphere.  The  distribution  functions,  or  scale  heights,  of  each  of  the  constituents  of  the 
heterosphere  are  found  b>  equating  the  pressure  gradient  of  the  atmosphere  with  the 
gravitational  force,  as  described  by  the  ideal  gas  law, 

P~(ty  <» 


where  p  =  gas  pressure,  k  =  Boltzman  constant,  m  =  molecular  mass,  T  =  temperature, 
and  p  =  density.  For  a  small  cross  sectional  area  of  thickness  dh 


therefore 


dP = -pgdh 

(2) 

dP  _  dP+dT 

P  p  T 

(3) 

^,=  *L(±+*Z) 

mg  [p  T  ) 

(4) 

By  assuming  isothermal  variation  and  defining  the  scale  height  to  be  H  =  — ,  a  simple 
differential  equation  is  obtained  described  by  the  following; 


&  =  ~dh  (5) 

P  H 

with  the  solution  being; 

p  =  p„e  "h  (6) 

Each  molecule  will  have  a  different  scale  height  depending  upon  its  mass.  This  gives  rise 
to  diffusive  equilibrium,  in  that  the  density  of  the  varying  constituents  will  decrease  at 
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certain  levels.  The  diffusion  process  takes  place  amongst  Ar,  O2,  N2,  O,  He,  and  H 
respectively  as  altitude  increases.  (Adler,  1992,  p.  1 1)  Figure  4  (  Akasofti,  1972,  p.  110) 
illustrates  the  densities  of  the  various  atmospheric  constituents  versus  altitude. 

7.  Time  Dependent  Variations 

Several  time  dependent  density  variations  are  present  in  the  earth's  atmosphere. 
The  density  of  the  earth's  atmosphere  varies  according  to  the  time  of  day,  day  of  the  year, 
and  which  year  it  happens  to  be  in  the  1 1  year  solar  cycle. 

a.  Diurnal  Variation 

The  atmospheric  density  variation  which  is  dependent  upon  the  time  of  day  is 
called  the  diurnal  variation.  The  maximum  density  of  the  upper  atmosphere  can  be  found 
at  approximately  1400  local  time,  with  the  minimum  around  0200  local  time.  The  density 
variation  becomes  more  pronounced  with  an  increase  in  altitude  as  can  be  seen  in  Figure  5. 
(Hess,  1968,  p.  99)  The  diurnal  variation  is  caused  by  the  alternate  heating  during  the  day 
and  cooling  during  the  night  of  the  upper  atmosphere.  Often  called  the  diurnal  bulge,  the 
density  variation  occurs  mainly  over  the  equator,  with  an  elongation  in  the  north-south 
direction  due  to  the  tilt  of  the  ecliptic.  The  peak  of  the  phenomena  occurs  at  the  latitude 
of  the  sub-solar  point.  (Jacchia,  1963,  p.  983)  The  temperature  change  in  the  upper 
atmosphere  during  this  diurnal  correction  parallels  that  of  the  density,  except  that  the 
temperature  change  lags  behind  the  density  change  by  approximately  two  hours.  This  lag 
seems  opposite  of  what  would  be  expected,  and  is  not  completely  understood. 

b.  27-day  Variations 

It  has  been  established  that  there  is  a  density  variation  that  is  caused  by  the 
27-day  solar  cycle.  This  was  shown  by  Jacchia,  who  determined  the  correlation  between 
actual  satellite  drag  measurements  and  the  solar  decimetric  flux,  or  10.7  cm  flux.  It  was 
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found  that  the  density  of  the  atmosphere  not  only  had  a  daily  variation,  but  also  a  monthly, 
or  27  day  cycle,  that  corresponded  to  the  27  day  rotation  cycle  of  the  sun.  (Hess,  1968, 

P  101) 

c.  Semi-Annual  Variations 

One  of  the  least  understood  upper  atmosphere  density  variations  is  the  semi¬ 
annual  variation.  This  density  variation  is  characterized  by  a  pronounced  minimum  during 
the  June-July  time  frame,  with  a  secondary  minimum  occurring  in  January.  The 
maximums  occur  during  the  September-October  time  frame,  with  a  lesser  secondary 
maximum  occurring  during  March-April.  Several  theories  exist  as  to  what  causes  the 
semi-annual  variation.  The  most  controversial,  is  that  the  variation  is  an  effect  of  the  solar 
wind.  Another  theory  is  that  the  variation  is  caused  by  the  convective  flows  from  the 
summer  pole  to  the  winter  pole.  The  flows  would  descend  at  the  winter  pole,  transporting 
heat  to  the  cooler  mesosphere  from  the  higher  altitudes.  (Hess,  1968,  p.  103) 

(L  Long  Term  Variations 

It  was  found  that  there  was  a  long-term  density  variation  associated  with  the 
1 1  year  solar  cycle.  Measurements  of  the  solar  flux,  or  FI  0.7,  taken  over  several  years 
provide  evidence  of  this  effect.  As  can  be  seen  in  Figure  6  (Larson,  1992,  p.  209),  the  sun 
has  an  1 1  year  cycle  with  maximum  and  minimum  values  of  F10.7.  The  F10.7 
measurement,  explained  in  detail  later,  is  a  measurement  of  the  EUV  radiation.  During 
periods  of  high  solar  activity,  the  solar  flux  is  high,  thus  causing  an  increase  in  atmospheric 
density  due  to  EUV  heating  and  the  resultant  expansion.  As  can  be  assumed  from  the 
figure,  the  density  of  the  upper  atmosphere  has  a  corresponding  eleven  year  minimum  and 
maximum. 
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B.  THERMOSPHERIC  PROCESSES 

The  upper  atmosphere  or  thermosphere  undergoes  a  change  in  composition  and 
density  due  to  several  external  inputs.  The  majority  of  the  change  is  caused  in  response  to 
the  activity  of  the  sun.  The  radiation  from  the  sun,  both  thermal  and  ultraviolet,  cause 
varying  rates  of  change  in  the  composition  and  hence  the  density  of  the  atmosphere.  The 
sun  also  causes  changes  in  density  due  to  the  solar  wind.  The  impingement  of  the  solar 
material  upon  the  earth's  magnetosphere  and  upper  atmosphere  causes  several  changes  in 
the  make-up  of  the  thermosphere  and  subsequently  the  overall  density  at  altitude.  The 
other  major  contributor  to  density  variation  of  the  earth's  atmosphere  is  geomagnetic 
activity  Geomagnetic  storms,  although  short-lived,  cause  a  significant  change  in  the 
atmospheric  composition  and  density. 

1.  Solar  EUV  Radiation  Effects 

The  first  source  of  density  change  that  will  be  addressed  is  the  Extreme 
Ultraviolet  Radiation  produced  by  the  sun.  The  sun's  EUV  radiation  is  the  main  cause  of 
density  variation  in  the  earth's  thermosphere.  The  solar  EUV  radiation  is  deposited  mainly 
at  low  latitudes  in  the  summer  hemisphere.  (Marcos,  1993,  p.  2)  The  circulation  and 
structure  of  the  thermosphere  at  the  low  and  middle  latitudes  are  controlled  by  the  heating 
caused  by  the  absorption  of  the  EUV  radiation  by  the  ozone  layer  in  the  lower 
thermosphere.  The  longer  wavelength  UV  and  visible  radiation  reach  the  lower 
atmosphere  and  hence  heat  the  earth's  surface.  (Bums,  1991,  p.  3)  Most  of  the  EUV 
radiation  reaching  the  earth's  atmosphere  is  absorbed  at  the  300  km  level  and  the  energy 
enters  the  atmosphere  through  photoionization.  The  energy  that  has  been  absorbed  by  the 
electrons  and  ions  is  passed  to  the  neutral  atmosphere  by  collisional  processes.  (Bums, 
1991,  p.3) 
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The  solar  EUV  radiation  also  imparts  a  momentum  to  the  neutral  gas  by  the 
creation  of  pressure  gradient  forces  that  drive  the  neutral  winds  from  the  day  to  night 
regions  and  from  the  summer  to  winter  hemisphere.  (Bums,  1991,  p.  3)  Variations  in  the 
strength  of  the  EUV  radiation  interacting  with  the  thermosphere  lead  to  changes  in  the 
composition  and  density  of  the  neutral  atmosphere.  During  the  period  of  solar  minimum, 
the  EUV  radiation  produced  by  the  sun  is  much  less  than  that  being  produced  during  solar 
maximum.  Therefore,  the  thermospheric  temperatures  and  neutral  densities  will  be  lov 
during  a  solar  minimum  than  during  solar  max.  This  fact  is  illustrated  in  Figure  7.  (Lar. 
1992,  p.  209)  It  can  be  seen  from  this  figure  that  the  variation  in  density  becomes  much 
greater  at  higher  altitudes  during  periods  of  high  solar  activity  versus  low  solar  activity. 

Due  to  the  fact  that  the  solar  EUV  radiation  is  absorbed  by  the  thermosphere,  it 
makes  it  difficult  to  obtain  a  measurement  of  the  EUV  impinging  on  the  atmosphere.  This 
problem  was  solved  by  Jacchia  in  the  early  sixties.  Jacchia  discovered  that  the  intensity 
found  at  10.7cm  closely  corresponded  to  the  amount  of  solar  activity  being  witnessed. 
Currently  the  accepted  measurement  of  solar  flux  is  the  F10.7  index.  As  can  be  seen  in 
Figures  6  and  7,  the  typical  value  for  FI 0.7  during  solar  minimum  is  75,  and  has  reached 
as  high  as  290  during  the  solar  maximum  period  of  1 958.  A  change  in  solar  activity  of  this 
proportion  would  mean  a  change  in  density  by  a  couple  orders  of  magnitude  at  the  altitude 
regime  of  interest. 

One  of  the  drawbacks  of  using  FI  0.7  as  a  gauge  of  EUV  radiation,  is  the  fact 
that  it  lies  at  the  other  end  of  the  spectrum  from  the  EUV  radiation,  and  is  not  a  direct 
measure  of  the  amount  of  EUV  radiation  reaching  the  thermosphere.  Figure  8  (Hess, 

1968,  p.  668)  illustrates  the  solar  spectrum.  It  can  be  seen  that  the  EUV  radiation  is  at  a 
frequency  on  the  left  hand  side  of  the  peak  spectrum,  whereas  the  solar  flux  measurement, 
F10.7,  is  on  the  right.  Because  the  F10.7  index  is  not  a  direct  measurement  of  the  amount 
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of  solar  EUV  radiation  entering  the  atmosphere,  accurate  density  calculations  based  on 
solar  activity  have  some  inherent  error.  Several  programs  have  been  initiated  in  order  to 
remedy  this  problem,  but  currently  the  FI  0.7  index  is  the  best  indicator  of  solar  flux 
available,  and  consequently  is  the  index  most  widely  used  in  atmospheric  modeling. 

2.  Solar  Wind 

Another  of  the  driving  forces  causing  composition  and  density  variation  in  the 
thermosphere  is  the  solar  wind.  The  solar  wind  consists  of  protons  and  electrons  flowing 
outward  from  the  sun's  corona.  Higher  density  plasma  streams  are  also  ejected  from  the 
sun  during  periods  of  flare  and  sunspot  activity.  (Fleagle,  1963,  p.  236)  The  solar  wind 
blows  the  interplanetary  magnetic  field  lines  across  the  polar  cap  in  a  direction  away  from 
the  sun.  (Bums,  1991,  p.  4)  This  in  turn  causes  a  potential  drop  across  the  earth's 
magnetic  polar  cap  as  the  interplanetary  magnetic  field  encounters  the  earth's  magnetic 
field.  Field-aligned  currents  flow  down  to  the  ionosphere,  closing  the  circuit,  and 
producing  an  ion-convection  pattern  at  high  latitudes.  The  ions  in  this  convection  pattern 
collide  with  the  neutral  particles,  driving  them  in  a  similar  convection  pattern.  (Bums, 
1991,  p.  4)  These  collisions  produce  heat,  which  in  turn  produces  a  rising  motion  around 
the  auroral  zone.  The  up  welling  and  the  convection-driven  neutral  winds  produce  a 
composition  and  density  change  which  spreads  from  the  high  latitudes  to  the  lower 
latitudes. 

The  convection  driven  neutral  winds  also  produce  another  significant  heat 
source.  Joule  heating  results  from  the  frictional  heating  of  the  plasma  as  it  is  dragged 
through  the  neutral  upper  atmosphere  by  the  auroral  electric  field  forces.  (Marcos,  1993, 
p.  3)  This  Joule  heating  is  a  substantial  heat  source,  but  becomes  even  more  prevalent 
during  periods  of  solar  flare  activity.  Flares  on  the  sun  cause  the  solar  wind  to  accelerate, 
driving  the  interplanetary  magnetic  field  faster  across  the  earth's  magnetic  field.  This 
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increases  the  cross-cap  potential  and  the  rate  of  panicle  precipitation,  and  can  also 
produce  a  magnetic  substorm.  (Bums,  1991,  p.  5) 

3.  Geomagnetic  Storms 

Another  of  the  major  contributors  of  composition  and  density  variation  are 
geomagnetic  storms.  Geomagnetic  storms  are  produced  by  the  interaction  of  the 
interplanetary  magnetic  field  with  the  earth's  magnetic  field.  When  solar  flare  activity  is 
occurring,  shock  waves  in  the  interplanetary  magnetic  field  are  driven  into  the  earth's 
magnetic  field  causing  rapid  transients  in  the  earth's  magnetic  field.  This  effect  is 
monitored  by  means  of  the  Ap  index,  or  geomagnetic  activity  index,  during  periods  of 
solar  flare  activity.  The  geomagnetic  storm  presents  itself  as  a  world-wide  transient 
variation  in  the  earth's  magnetic  field.  The  onset,  or  sudden  commencement,  of  a 
magnetic  storm  is  characterized  by  a  rapid  increase  in  the  Ap  index.  Approximately  20 
minutes  later,  the  temperature  and  density  of  the  auroral  zones  begins  to  increase. 

One  of  the  manifestations  of  geomagnetic  storms  is  the  generation  of  waves 
that  propagate  from  the  auroral  regions  to  the  lower  latitudes.  These  waves  take 
approximately  eight  hours  to  reach  the  lower  latitudes.  (Alder,  1992,  p.  18)  A  typical 
magnetic  storm,  illustrated  in  Figure  9  (Akasofu,  1972,  p.  557),  lasts  from  24  to  48  hours 
or  longer.  The  effects  of  the  geomagnetic  storm  on  the  density,  last  even  longer  and  are 
quite  pronounced  at  the  higher  altitudes.  This  is  illustrated  in  Figure  10.  (Ratcliffe,  1972, 
p.  35)  It  can  be  seen  from  Figure  10,  that  during  periods  of  geomagnetic  storms,  the 
density  at  the  altitude  regime  of  interest  increases  several  fold.  This  fact  coupled  with 
intense  bombardment  of  EUV  radiation,  increases  the  density  by  several  orders  of 
magnitude. 
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4.  Gravity  Waves,  Planetary  Waves,  and  Tides 

The  final  source  of  composition  and  density  variation  to  be  discussed  are 
propagating  tides  and  gravity  waves.  Atmospheric  tides  are  caused  primarily  by  the 
absorption  of  ultra-violet  radiation  by  the  ozone  in  the  stratosphere,  while  gravity  waves 
are  caused  by  a  variety  of  mechanisms  which  occur  in  the  troposphere.  A  couple  of  these 
mechanisms  are  the  shears  associated  with  cold  fronts  and  winds  blowing  over  mountains. 
(Bums,  1991,  p.  3)  At  lower  altitudes,  the  semi-diurnal  tide  is  the  major  contributor  to 
density  variation.  This  effect  becomes  less  apparent  at  higher  altitudes  due  to  the  fact  that 
the  semi-diurnal  tide  is  dissipated  by  viscosity  and  ion  drag.  At  higher  altitudes,  the 
diurnal  tide  becomes  the  driver  of  density  variation  and  is  caused  by  the  absorption  of  the 
EUV  radiation  in  the  thermosphere.  Overall,  these  waves  and  tides  propagate  up  from  the 
lower  altitudes  affecting  the  composition  and  density  of  the  upper  altitudes. 

As  a  consequence  of  the  above  mentioned  thermospheric  processes,  it  is  known 
that  the  temperature  and  hence  the  density  of  the  upper  atmosphere  vary  with  the 
following  parameters: 

•  solar  flux  (solar  cycle  and  daily  component) 

•  geomagnetic  activity 

•  local  time 

•  day  of  the  year 

•  latitude 

•  longitude 

•  wave  structures 
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C.  PROGRAM  DEVELOPMENT 

The  propagation  program  created  for  the  atmospheric  model  comparison  is  a 
FORTRAN  based  program  which  uses  a  Runge-Kutta  integrator.  The  propagator  uses  a 
form  of  Cowell's  Method  of  orbital  prediction.  The  program  uses  the  position  and 
velocity  vectors  in  Cartesian  coordinates  as  its  inputs,  and  integrates  over  a  designated 
time  frame  to  produce  the  position  and  velocity  vectors,  also  in  Cartesian  coordinates. 
Cowell's  Method  of  orbital  prediction  has  been  found  to  be  inaccurate  due  to  the  build-up 
of  round  off  error.  A  more  accurate  method  would  have  been  to  convert  the  Cartesian 
coordinates  into  normalized  spherical  coordinates  in  order  to  minimize  the  integration 
error  build-up.  Further  research  into  this  conversion  is  being  conducted  in  order  to  obtain 
greater  accuracy.  For  the  purpose  of  the  comparison,  however,  the  accuracy  of  the  result 
was  not  in  question,  but  the  accuracy  of  the  atmospheric  model,  and  its  relative  ease  of 
use. 

The  propagation  program  consists  of  a  series  of  subroutines  which  calculate  the 
perturbing  forces  acting  on  the  satellite.  These  forces,  in  the  form  of  accelerations,  are 
applied  to  the  satellite's  motion  resulting  in  a  position  and  velocity  vector  reflecting  the 
result  of  the  perturbing  forces.  Figure  1 1  and  Table  l(Milani,  1987,  p.  14-15)  illustrate 
the  various  perturbing  accelerations  and  their  relative  magnitude  as  a  function  of  altitude. 
In  order  to  obtain  an  accurate  reflection  of  orbital  motion  at  the  altitude  band  in  question, 
it  was  decided  that  in  addition  to  the  drag  force,  the  earth's  geopotential  and  third  body 
attractions  would  also  be  included  as  perturbing  forces. 

1.  Earth's  Geopotential 

As  can  be  seen  from  Figure  1 1,  the  main  perturbing  force  encountered  by  low 
earth  orbiting  satellites  is  the  earth's  gravity  field.  If  the  earth  was  a  perfectly  round  and 
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smooth  planet,  then  the  forces  of  gravity  could  easily  be  modeled  by  the  inverse  square 
gravity  law  shown  in  equation  7. 

u 

S  =  “  3  r  (7) 

r 

However,  since  the  earth  is  not  perfectly  spherical,  and  not  smooth,  the  gravity  field  must 
be  derived  by  obtaining  the  solution  to  Laplace's  equation: 

v2V  =  0  (8) 


where 


_  r  dm 

V  =  G  I  —  (9) 

volume 

8  being  the  distance  from  the  satellite  to  an  incremental  mass,  dm,  inside  the  earth,  and  G 
is  the  factor  of  proportionality  in  Newton's  Law  of  Gravitation.  Laplace  then  derived  the 
basic  equation  for  the  geopotential  shown  below: 


By  converting  to  spherical  coordinates  and  applying  Rodriguez’  formula,  the  earth's 
geopotential  takes  the  form  of  equation  1 1 . 


V  =  —  1  +  XXf-)  Pnm{sm<p'){Cn„cosm\  +  sinm\)  (11) 

In  this  equation,  V  is  the  gravitational  potential  function,  GM  is  the  earth's  gravitational 
constant,  r  is  the  radius  vector  from  the  earth's  center  of  mass,  a  is  the  semi-major  axis  of 
the  model  ellipsoid,  n  and  m  designate  the  degree  and  order  of  the  coefficients,  <p  is  the 
geocentric  latitude,  X  is  the  geocentric  longitude,  C„m  and  Snm  are  the  harmonic 
coefficients,  and  Pnm  are  the  associated  Legendre  functions.  (Ross,  1993,  p.5-3)  Several 
Earth  gravitational  models  are  currently  in  use  in  both  the  civilian  and  DoD  communities. 
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The  DoD  currently  uses  the  WGS-84  geopotential  model  as  its  gravitational  model.  The 
WGS-84  Earth  gravitational  model  (EGM)  contains  a  180  by  180  matrix  of  zonal  and 
tesseral  harmonic  coefficients.  The  WGS-84  EGM  is  currently  being  used  for  several 
DoD  satellite  programs  including  the  GPS  satellite  constellation,  and  is  considered  to  be 
an  extremely  accurate  representation  of  the  earth's  geopotential.  For  most  cases,  the  full 
matrix  is  not  needed,  and  the  matrix  is  truncated  down  to  a  41  by  41  matrix.  The 
truncated  matrix  provides  an  ample  representation. 

Currently,  the  propagation  program  contains  the  first  8  sets  of  zonal  and  tesseral 
coefficients.  Future  iterations  of  the  propagation  program  will  contain  the  complete  41  by 
41  matrix.  At  the  moment,  attempts  at  configuring  the  propagation  program  in  the  WGS- 
84  coordinate  system  have  failed.  The  first  eight  sets  of  coefficients  vary  only  slightly 
from  the  first  eight  zonal  terms  of  the  more  basic  geopotential  models  and  can  be  used 
without  incurring  significant  errors  during  propagation. 

2.  Drag 

For  satellites  in  low  earth  orbit,  drag  is  the  other  major  perturbing  force  that 
must  be  modeled.  Drag  affects  all  satellites  in  all  altitude  regimes,  but  the  affect  is 
considered  to  be  insignificant  at  altitudes  greater  than  1000  km.  The  following  equation 
represents  the  atmospheric  drag  acceleration  which  is  placed  on  an  orbiting  body. 

a°=\C°ZpV‘  (12) 

i  m 

In  this  equation,  aD  represents  the  drag  acceleration  imparted  upon  an  orbiting  satellite, 

CD  is  the  coefficient  of  drag,  A  is  the  cross  sectional  area  of  the  satellite  perpendicular  to 
the  direction  of  motion,  m  is  the  satellite  mass,  V  is  the  satellite  velocity  and  p  is  the  local 
atmospheric  density.  The  velocity  used  in  this  equation  is  computed  by  combining  the 
geocentric  velocity  of  the  satellite  with  the  contribution  due  to  earth  rotation  and  the  wind 
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velocity  at  altitude.  It  is  important  to  note  that  the  wind  component  of  the  velocity  cannot 
be  ignored,  and  can  be  quite  significant  at  altitudes  in  excess  of  600  km.  The  neutral  wind 
is  a  consequence  of  the  activity  caused  by  a  geomagnetic  disturbance.  The  change  in  drag 
can  be  5%  for  every  200m/s  of  wind  velocity.  During  a  period  of  intense  geomagnetic 
activity,  the  neutral  winds  can  reach  velocities  in  excess  of  1  km/s  which  is  equivalent  to  a 
25%  change  in  density.  (Marcos,  1993,  p.  6)  In  order  to  model  this  behavior,  several 
atmospheric  wind  models  have  been  developed.  Models  of  the  neutral  winds  begin  with 
the  efforts  of  Sissenwine  in  the  late  sixties,  to  the  Horizontal  Wind  Model  90  (HWM  90). 
Currently  these  wind  models  are  not  incorporated  into  the  empirical  atmospheric  drag 
models,  but  are  incorporated  in  the  general  circulation  models. 

The  coefficient  of  drag,  CD,  is  a  difficult  quantity  to  obtain.  In  order  to 
determine  the  coefficient  of  drag,  it  must  be  determined  whether  the  satellite  is  in  a 
continuum  flow,  or  a  free  molecular  flow.  This  is  done  by  determining  the  Knudsen 
number.  The  Knudsen  number  is  the  ratio  between  a  typical  dimension  of  the  satellite  and 
the  average  mean-free-path  of  the  molecules  found  in  the  local  atmosphere.  When  the 
Knudsen  number  is  much  less  than  one,  the  satellite  is  considered  to  be  in  a  continuum 
flow.  The  satellite  is  considered  to  be  in  a  free  molecular  flow  when  the  Knudsen  number 
is  greater  than  one.  Since  the  atmospheric  density  is  so  low  at  orbital  altitudes,  satellites  in 
the  upper  atmosphere  are  characterized  by  free-molecular-flow  aerodynamics.  (Ross, 

1994,  p.  16) 

The  collision  of  the  atmospheric  particles  with  the  spacecraft  produce  the 
atmospheric  drag  force.  The  collisions  can  be  classified  under  three  categories:  (1)  elastic 
or  specular  reflection,  (2)  diffuse  reflection  and  (3)  absorption  and  subsequent  emission. 

In  elastic  or  specular  reflection,  the  molecule  collides  with  the  satellite  and  is  reflected 
away  without  transferring  energy  to  the  satellite.  In  diffuse  reflection,  the  molecule 
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collides  with  the  satellite  and  transfers  a  portion  of  its  energy.  Energy  is  also  transferred 
to  the  satellite  when  molecules  are  absorbed  on  impact  and  later  emitted.  Since  the 
satellite  is  considered  to  be  in  a  free-molecular-flow,  the  coefficient  of  drag  is  obtained 
from  a  statistical  mechanical  calculation.  (Ross,  1994,  p.  16)  Depending  on  the  size  and 
makeup  of  the  satellite,  the  coefficient  of  drag  can  vary  from  2.0  to  6.0.  Table  2  (Larson, 
1992,  p.  207)  illustrates  the  variation  of  the  drag  coefficient  for  various  orbiting  platforms. 
If  the  coefficient  of  drag  cannot  be  determined,  a  value  of  2.2  is  used.  For  the  purpose  of 
this  thesis,  the  coefficient  of  drag  will  remain  constant  at  a  value  of  2.2. 

In  equation  12,  the  coefficient  of  drag,  the  cross  sectional  area,  and  the  mass  of 
the  satellite  all  make  up  what  is  called  the  B-factor. 

B  =  ^-  (13) 

m 

Inverting  equation  13  will  result  in  what  is  more  widely  known  as  the  Ballistic  coefficient. 
Current  practice  is  to  guess  what  the  correct  B-factor  is  for  that  given  day  and  adjust  the 
B-factor  until  the  correct  position  and  velocity  is  achieved.  This  does  not  seem  to  be 
orbital  prediction,  but  rather  orbital  correction.  If  the  size  and  shape  of  the  satellite  are 
known,  as  well  as  the  mass,  then  the  B-factor  should  only  vary  with  varying  coefficients 
of  drag.  However,  most  of  the  current  orbital  prediction  schemes  use  the  B-factor  as  a 
catch-all  for  any  other  errors  in  the  modeling  program.  Hence  the  value  of  the  B-factor  is 
nowhere  near  the  actual  value.  In  this  comparison,  it  was  decided  to  keep  the  B-factor  at 
a  constant  value  in  order  to  obtain  a  more  accurate  comparison  between  atmospheric 
models.  It  must  be  noted  at  this  point  that  keeping  the  B-factor  at  its  actual  value  is  a 
major  change,  vice  using  the  B-factor  as  an  error  catch-all. 

The  density  for  the  drag  acceleration  calculations  is  of  course  derived  from  the 
atmospheric  models.  In  the  calculation  of  the  drag  acceleration,  the  density  is  the  most 
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difficult  to  obtain,  and  frequently  the  one  parameter  with  the  greatest  error.  As  a  rule,  the 
most  current  atmospheric  models  have  a  1 5  to  20%  inaccuracy  rate,  with  this  inaccuracy 
increasing  as  altitude  increases.  Due  to  this  inaccuracy,  the  propagation  scheme  is  only  as 
accurate  as  the  atmospheric  model  being  used. 

3.  Third  Body  Attractions 

The  last  of  the  perturbing  forces  currently  included  in  the  propagation  program  is 
that  of  third  body  attractions.  In  the  case  of  the  satellite  program  in  question,  the 
perturbing  bodies  are  the  sun  and  moon.  As  can  be  seen  in  Figure  1 1 ,  both  the  sun  and  the 
moon  contribute  some  small  portion  of  disturbance  force  to  a  medium  altitude  satellite.  In 
practice,  this  disturbance  force  is  usually  ignored  for  the  low  earth  orbiting  platforms,  but 
in  the  altitude  band  in  question  (600  -  800km),  these  disturbance  forces  must  not  be 
ignored  if  a  truly  accurate  representation  of  orbital  motion  is  to  be  modeled.  The  equation 
used  to  find  the  perturbing  acceleration  due  to  the  moon  is  described  by  equation  14 
below. 


In  this  equation,  r  is  the  radius  from  the  earth  to  the  satellite,  rm  is  the  radius  fiuir.  the 
moon  to  the  satellite,  r ^  is  the  radius  from  the  earth  to  the  moon,  and  is  the 


gravitational  parameter  of  the  moon.  (Bate,  1971,  p.  389) 

Currently  the  propagation  scheme  does  not  contain  several  perturbations  that  are 
important  for  accuracy  purposes.  At  the  moment,  the  sun’s  radiation  pressure  is  ignored 
as  well  as  the  earth's  albedo  effect.  Also,  relativistic  effects  are  ignored,  which  must 
eventually  be  incorporated  in  order  to  improve  accuracy. 
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D.  ATMOSPHERIC  MODEL  DEVELOPMENT 

As  mentioned  previously,  atmospheric  models  can  be  divided  into  two  categories, 
the  empirical  models,  and  the  general  circulation  models.  The  general  circulation  models 
are  dynamic  representations  of  the  earth's  atmosphere,  and  require  extensive 
computational  time  in  order  to  model  the  atmosphere.  The  general  circulation  models  that 
have  been  recently  created  require  Cray  computers  to  run  simulations.  Efforts  are 
currently  under  way  to  convert  these  models  to  a  more  user  friendly  format,  so  that  they 
may  be  used  on  smaller  computers.  The  empirical  models,  however,  are  readily  available, 
and  take  little  computational  time  per  simulation.  One  drawback  is  that  the  accuracy  of 
the  models  is  quite  a  bit  less  than  that  of  the  general  circulation  models. 

The  history  of  the  empirical  earth  atmospheric  models  dates  back  to  the  efforts  of  L. 
G.  Jacchia  in  the  late  fifties.  Once  the  early  satellites  such  as  Sputnik  and  Pioneer  had 
been  launched,  immediate  drag  analysis  was  performed,  and  atmospheric  models 
developed  based  on  this  analysis.  The  early  models  were  crude,  and  only  represented  the 
regions  where  the  satellites  were  orbiting.  Little  was  understood  of  the  variability  of  the 
environment  and  the  density  fluctuations  being  encountered.  As  more  satellites,  and  a 
greater  understanding  of  the  sun's  interaction  with  the  earth's  atmosphere  was  obtained, 
the  accuracy  of  the  atmospheric  models  increased.  Figure  12  illustrates  the  developmental 
history  of  the  various  earth  atmospheric  drag  models.  (Marcos,  1993,  p.  20)  The 
mathematical  development  of  the  Lockheed-Densel  or  Jacchia  60  model,  the  Jacchia  71, 
and  the  MSIS  86  earth  atmospheric  models  are  described  below. 
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1.  Lockheed  Densel 

The  Lockheed  Densel,  or  Jacchia  60  earth  atmospheric  model  was  the  first 
model  to  be  implemented  into  the  prediction  scheme.  As  can  be  seen  from  Figure  12,  it 
was  one  of  the  earliest  atmospheric  models  contrived,  and  hence  its  accuracy  is  in  great 
question.  The  Lockheed  Densel  model  is  actually  a  combination  of  two  atmospheric 
models;  the  ARDC  1959  density  model  for  low  altitudes  (h  <  76  nautical  miles),  and  the 
Jacchia  60  density  formulation  for  h  >  76  nautical  miles.  For  the  density  below  76 
nautical  miles,  the  density  is  obtained  from  the  ARDC  1959  model  which  contains  a  table 
of  density  values  as  a  function  of  altitude.  The  discussion  on  obtaining  the  density  below 
76  nm  will  not  be  discussed  in  this  paper.  When  the  density  at  an  altitude  above  76  nm  is 
desired,  the  Jacchia  60  formulation  is  used.  The  first  requirement  is  to  define  the  unit 
vector  to  the  diurnal  bulge  using  the  solar  position  and  the  bulge  lag  angle.  In  order  to 
accomplish  this  the  solar  longitude  must  be  determined  by  the  following  equation, 

''•=(2%5.25)-I-4l+00335si"(2%5.25)  <l5> 

where  d  represents  the  number  of  days  since  December  31,  1957.  The  unit  vector  to  the 
sun  is  obtained  from  the  following  series  of  equations. 

t  =  cos^  m.  =  sin  A.  cose  n,  =  sin  A-  sin  e  (16) 

where  e  is  the  obliquity  of  the  ecliptic.  (Lockheed,  1992,  p.  B-100) 

The  unit  vector  to  the  diurnal  bulge  is  then  calculated  in  true  of  date  coordinates 
by  the  following  matrix 
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UB  =  C 


&cos0-m.sin0 
m,cos0+Lsin  6 

th 


(17) 


where  C  =  J2000.0  to  true  to  date  transformation  matrix,  and  8  =  0.55  radians  which 
equals  the  bulge  lag  angle.  Two  options  are  available  to  the  user  when  using  the  solar  flux 
value  of  FI  0.7.  The  user  may  either  specify  a  value,  or  one  is  calculated  using  the 
following  formula, 

F..,  =  15+8cos(2®^020)  (18) 

where  once  gain  d  is  the  number  of  days  since  December  31,  1957.  This  equation  allows 
an  approximation  of  the  F10.7  value  on  any  given  day  over  the  1 1-year  solar  cycle  period. 
(Lockheed,  1992,  p.  B-103)  The  Jacchia  60  model  divides  the  atmosphere  into  a  series  of 
three  bands  for  density  calculation.  The  first  band  is  from  76-  108  nautical  miles.  The 
equation  used  to  calculate  the  density  in  this  altitude  band  is  given  by 


X[10+[l2^)10+COS^3 

p  =  7.18 

(p)™  =  density  from  ARDC  1959  model  at  76  nm 
h  =  altitude  in  nm 
FW1  -  solar  flux  measurement 
R*U, 


(19) 


cos<p  = 


R 


R  =  SV  true  of  date  position  vector 


(Jt  =diumal  bulge  vector  (equation  17). 
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The  second  altitude  band  ranges  from  108  -  378  nautical  miles,  and  the  density 
calculation  is  given  by  equation  20; 

p(h)  =  p.{h) 0. 85Fio  7[l.  0 + 0. 023  75(e00102*  -  1 . 9)(  1. 0 + cos<p)3  ]  (20) 

where  p.(/i)  =  k  exp|(-6  -  ab+  6. 363e‘0004®,, )  log.  1  oj 

a  -  0. 00368,  b  =  15.738,  and  k  is  the  conversion  factor  from  slugs  to  kg/km3  The  third 
and  final  altitude  band  that  is  calculated  in  the  Jacchia  60  model  lies  between  378  and 
1000  nautical  miles.  Equation  21  is  used  to  calculate  the  density  in  this  region. 

p{h)  =  (°  00504  125(1.0  + cos  ^)J(A3  -  6.0x  10‘)  +  6.0x  106J/r  <21> 

In  the  Jacchia  60  model  it  is  assumed  that  the  density  is  zero  when  the  altitude  is  above 
1000  nautical  miles. 


a.  Mathematical  Basis 

In  order  to  provide  a  mathematical  background  for  the  above  density  equations, 
the  partial  derivative  equations  are  illustrated  below.  The  partial  derivative  of  the 
calculated  density  with  respect  to  position  is  described  by  equation  22; 

dM  _  k  dfjh)  ^  dh  t  dp(h)  *  d(p  (22) 

dR  dh  dR  dtp  *  dR 

where  k\  is  the  conversion  factor  between  nautical  miles  and  kilometers.  The  partial 
derivative  of  the  altitude  with  respect  to  R  is  best  estimated  by  the  following  equation. 


R  =  position  vector  magnitude 


R.  =  earth  equatorial  radius 
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e  —  earth  eccentricity 
The  differentiation  yields  equation  24; 


<P'  =  geocentric  latitude 


dh  RT  n 
dR  R 


y/\  -e2e2  cos (/>' 
(l-e:  cos:  0')2 


dcosQ' 


dR 


(24) 


where 


X,  Y,  and  Z  represent  the  geocentric  coordinates  of  the  orbiting  object.  The  partial  of  <p 
with  respect  to  R  is  described  by  the  following  relationship. 


d(p  _  1 

dR  sin  <p 

The  partial  derivative  of  p(h)  with  respect  to  h  and  the  partial  of  p(h)  with 
respect  to  <p  are  different  for  each  of  the  altitude  bands.  In  the  altitude  band  between  76 
and  108  nautical  miles,  the  partial  derivative  of  p(h )  with  respect  to  h  and  (p  are  given  by 
the  following  equations: 


rR»U.-=r  U. 

- - — K - 

Ri  R 


(25) 


M/t) 

dh 


■Ppu-^C 

hy  32 


1  —  1. 1 33/**io  7 


+T224[l+“S,>]i  (2S) 

where  p(h )  is  the  density  found  from  equation  19,  p  =  7. 18,  and  the  variables  A,  B,  and  C 
follow  (Lockheed,  192,  p.  B-105): 


c  =  [10+f7^^1  +  cos<1P^ 

l  1224  J 


^^  =  -3AB(^-^\l  +  cos(py  sin  9  (27) 

d<p  V  1224  J 

For  the  altitude  band  between  108  and  378  nautical  miles,  these  two  equations 
take  on  this  form; 

=  -p(/»)[a  +  0. 0305424  exp(-0.0048A)]fw!0 

+0.0002059 1 25p.(h)F\o  ?[exp(0. 0 1 02 /»)]( 1  +  cos <p)3  (28) 

=  -0.0605625p.(A)Fio  7[exp(0.0102/i)-  1.9](l  +  cos<p)2  sin  <p  (29) 
d(p 

In  the  highest  altitude  band  between  378  and  1000  nautical  miles,  the  equation 
take  on  the  form  (Lockheed,  1992,  p.  B-106) 

^W=_?pW+2il21*2filZ.ti+cos®])t  (30) 

dh  h  h6  1 

jpW  =_ 0. 00 1 89f,._.r(l  +  sin ^’-6.0x10°)]*  (31) 

d<p  "  *■ 

These  equations  can  be  used  to  implement  an  atmosphere  in  an  orbital 
prediction  program,  and  in  fact  are  currently  used  for  that  purpose.  The  Jacchia  60 
density  model  was  one  of  the  first  empirical  models  on  the  market.  With  the  launching  of 
numerous  satellites  throughout  the  recent  years,  Jacchia  et.  al.  have  been  able  to  expand 
on  their  empirical  model.  As  can  be  seen  from  figure  1 2,  several  Jacchia  models  have  been 
developed,  each  model  building  on  its  predecessor. 


2.  Jacchia-Roberts  71 

The  next  model  to  be  discussed  is  the  Jacchia-Roberts  71  atmospheric  model. 
The  Jacchia-Roberts  71  atmospheric  model  takes  into  account  both  the  solar  and  the 
geomagnetic  activity  during  the  time  in  question.  L.  G.  Jacchia  defined  two  empirical 
profiles  to  represent  temperature  as  a  function  of  altitude  and  exospheric  temperature. 
One  profile  was  defined  between  the  region  of  90  to  125  kilometers,  and  the  other  above 
125  kilometers.  Jacchia  then  used  these  temperature  functions  in  the  appropriate 
thermodynamic  differential  equations  to  obtain  density  as  a  function  of  altitude  and 
exospheric  temperature.  (Lockheed,  1992,  p.  B-81)  The  Jacchia  model  as  it  stood  was 
very  cumbersome  and  required  a  great  amount  of  data  storage  to  hold  the  data  required, 
so  C.  E.  Roberts  provided  a  method  for  evaluating  the  Jacchia  model  analytically.  This 
lead  to  a  faster  and  more  manageable  model.  The  following  are  the  equations  used  in 
determining  the  atmospheric  density  as  derived  in  the  Jacchia-Roberts  71  atmospheric 
model. 

The  first  step  in  the  model  is  to  calculate  the  nighttime  minimum  global 
exospheric  temperature  for  zero  geomagnetic  activity, 

7c  =  379° +3°. 24Fio  ?  + 1°. 3[Fio.  i - F\o  7]  (32) 

where  F 10  7=81  day  running  average  of  the  FI  0.7  centered  on  the  day  in  question.  F 10  7 
=  solar  flux  measurement  as  obtained  from  the  solar  observatory  at  Ottawa,  Canada. 
(Jacchia,  1970,  p.  16) 

The  value  of  the  nighttime  minimum  exospheric  temperature  is  then  used  in 
calculating  the  uncorrected  exospheric  temperature  as  follows, 
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where 


T,=  £  jl  +  0.3  sin"  0+ (cos"  n- sin"  0) cos30  *  |  (33) 


=  e  =  -\<p+S\  r=//-37°.0+6°.0sin(//+43o.0) , 

2  2 


-K< t< K 

&  =  the  sun's  declination 

0 =the  geodetic  latitude  of  the  satellite  in  true  of  date  coordinates  ( includes 
earth  flattening) 


H  =  1 80° .  Q< 


{SxX2-S2X}) 

SiX'+S.X, 

4s1x2-1s;xi|C0S 

(34) 


The  X  variables  are  the  components  of  the  unit  position  vector  of  the  satellite  in  true  of 
date  (TOD)  coordinates,  and  the  S  variables  are  the  components  of  the  unit  vector  to  the 
sun  in  TOD  coordinates.  (Lockheed,  1992,  p.  B-83) 

It  has  been  found  through  analyzing  the  effect  of  geomagnetic  activity  on  the 
atmosphere,  that  there  is  a  lag  of  approximately  6  to  7  hours  from  a  detection  of  a  density 
change  from  the  actual  geomagnetic  disturbance.  In  order  to  account  for  this  lag,  the 
value  of  Kp  is  obtained  for  a  period  of  6.7  hours  prior  to  the  integration  time  in  question. 
It  must  be  noted  that  the  Kp  value  only  exists  at  a  three  hour  resolution.  At  this  point  the 
correction  for  the  exospheric  temperature  is  calculated  using  the  following  formulas, 

AT-  =  28° .  0  Kp  +  0° .  03eKp  (Z>  200  kilometers) 

AT-  =  14° .OKp  +  0° .02eKp  (Z<  200  kilometers)  (35) 

The  corrected  exospheric  temperature  is  then 
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T-  —  T\  +  A7’- 

and  the  inflection  point  temperature  is  given  by  the  following  formula  (Jacchia,  1970,  p. 
21) 

f,  =  37 1 0 . 6673  +  0.51 8806  T--  294  0 . 3 505e^°°21622 T~  (36) 

These  values  for  the  temperatures  and  the  satellite  altitude  are  used  in  the  calculation 
presented  by  Roberts  for  the  atmospheric  density.  However,  a  number  of  corrections 
must  be  applied  due  to  several  atmospheric  processes  presented  by  Jacchia.  These 
corrections  will  be  presented  before  continuing  with  the  Roberts  calculations. 

One  of  the  corrections  deals  with  the  direct  effect  of  geomagnetic  activity  on  the 
density  below  200  kilometers.  This  correction  is  calculated  by  the  following  relation, 
(Alogp)o  =  0.012  Kp  + 1.2  x  10‘V'  (37) 

The  next  correction  deals  with  the  semi-annual  density  variation.  This  correction 
is  calculated  from  the  following  relations,  where  Z  is  the  altitude  in  kilometers. 

(Alog p)sa  ~  f{Z)g{t)  (38) 

where 

/(Z)  =  (5.876  x  10-7Z233’  +  0.06328)e“° 0028682  (39) 

g(/)  =  0.02835+[0.3817+0.I7829sin(2^iw+4.137)]sin(4n-RJ +4.259)  (40) 


zsa  =  0+  0. 09544 j|j  +  |sin(2a^  +  6.035) 


1.65 


(41) 


JD\m 

365.2422 


(42) 


where  JD»*  is  the  number  of  Julian  days  since  January  1,  1958.  (Lockheed,  1992,  p.  B 
84) 
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The  next  correction  deals  with  the  seasonal  latitudinal  variations  in  the  lower 
thermosphere.  Equation  43  represents  the  general  density  variation,  whereas  equation  44 
represents  the  correction  for  Helium  specifically. 

(Alogp)ir  =0.014(Z  — 90)e^°°n(Z  901  ^  sin(2;r0  + 1.72)  sin  p|  sin  p|  (43) 

—  f  sin  -0.35355  j  (44) 

e[  (4  2|&|J  J 

where  e  is  the  obliquity  of  the  ecliptic.  (Lockheed,  1992,  p.  B-84) 

Below  125  kilometers,  Roberts  uses  the  same  temperature  profile  as  Jacchia, 

7 tz)  =  r,+^r£c,r  (45) 

35  naO 

where 

d\=Tx-  To,  To  =  183°. 0/C  and  the  coefficients  follow 

Co  =  -89284375.0  Cz  -  -52687. 5km~2  C4  =  -0.Skm~i 

C.  =  3542400. 0km~'  Cz  =  340.5AwT3  (46) 

where  Tx  is  the  inflection  point  temperature  at  125  kilometers  calculated  by  equation  36. 
Roberts  then  substituted  the  temperature  profile  obtained  from  equation  45  into  the 
barometric  differential  equation  and  integrated  by  partial  fractions  to  obtain  the  following 
expression.  Equation  47  represents  the  density  found  in  the  altitude  band  between  90  and 
100  kilometers.  (Lockheed,  1992,  p.  B-85) 

p'(z)=^(w)F,‘exp(4F!)  t47) 

where  the  "o"  indicates  the  conditions  at  90  kilometers.  The  mean  molecular  weight  is 
calculated  by  the  following  equation, 


(Alogp)tf«  =  0.65 
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(48) 


where 


M{z)  =  ^AnZn 

n»0 


Ao  =  -435093.363387  As  =  1 1.043387545Aw'3 

Ai  =  28275. 5646391^7//"'  Aa  =  -0.08958790995*///^ 

A  2  =  -765. 33466 1  08*///  "2  /4»  =  0.00038737586*///"5 

/is  =  -0. 000000697444 knt^ 

These  constants  give  a  value  of  the  mean  molecular  weight  at  90  kilometers  of  28.82678 
which  is  close  to  the  sea  level  mean  molecular  weight  of  28.960.  (Lockheed,  1992,  p.  B- 
86)  The  density  of  the  lower  limit  is  assumed  to  be  constant  at  p.  =  3.46  x  10"9  gm/cm 3 
The  constant  k  in  equation  47  is  evaluated  by  the  following  formula 


3  s4gR;- 

RdiCA 


(49) 


where  g  is  assumed  to  be  9.80665m/sec2  ,  the  acceleration  due  to  gravity  at  sea  level. 
The  functions  Fx  and  Fz  in  equation  47  are  determined  from  equations  50  and  5 1 . 


Z+Rq 

90+/? 


Z-r, 


W 


90 -r, 


i  J 


Z-r , 


90 -r. 


2/ 


z2-ixz  +  x2+y 2 

8100-180Y  +  Y2+T2 


(50) 


K=(Z-  90) 


A,  + 


(Z+/?J(90+/?J 


+— tan"1 


r(z-9o) 


V2+(Z-X)(90-X) 


(51) 


The  variables  r,  and  r2  are  the  two  real  roots  and  X  and  Y  are  the  real  and 
imaginary  parts  (Y>0),  respectively,  of  the  complex  conjugate  roots  of  the  following 
quadratic. 


with  the  following  coefficients 


/>(z)=Xc;,z' 

m*Q 


(52) 
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l<n<4 


3  54  T  C 

and 

°  CA  C4 


c 

n 

"~C 


for  values  of  C„  used  in  equation  45.  The  parameters  p,  in  equations  50  and  5 1  are 
evaluated  by  the  following  relations  (Lockheed,  1 992,  p.  B-86) 


S(r,)  _  _-S(r,)  _  _S(-RJ 

A=C/(r,)  A '  t/(rj  K 

p.  =  {s0  -r,r,Rl[B,  +  B,( 2X+r,  +r2 -  Rj]+Bir,)p2}/ X' 

-{r,r,AR.{ X‘  +  >")+ "'('i  )ft  +V, (*.!  -X'-r1  )ft  }  /  A" 

ft  +B,(2JT  +r,  +r.  -R.)-p, -2(X  +  R„)p,  ~(r,  +  Rjp, -(r:+Rjp, 
pi=Bi-2p,-pJ-p2 

where 

X*  =  -2r/2Ra(R2a  +2  XRa  +  X2+Y2) 

V  =  (*,  +r,)to  +/;)to2  +2XRa  +  X2+Y2) 

U(r,)  =  (r,+Raf{r2 -2Xr, +X2  -h F2 )(r,  -r2) 

(W  1  Y1  ^ 

R.+- — — 

r<  J 

and 


B,  =  a.+fin-^—  <»  =  0,1 . 5) 

S(Z)  =  £fl„Z" 

n*0 

with  the  coefficients 

a0  =3144902516.672729  j30  =  -52864482.17910969 

a,  =  -123774885.4832917  /3,  =  -16632.50847336828 
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a,  =  1816141.096520398  /3,  =  -1.308252378125 

a3  =  -1 1403.31079489267  03  =  0.0 

a4  =24.36498612105595  04=O.O 

a,  =  0.008957502869707995  /3,  =  0.0 

Equation  47  is  a  valid  equation  below  1 00  kilometers  where  mixing  is  assumed 
to  be  predominant.  However,  above  100  kilometers,  diffusive  equilibrium  is  assumed,  and 
it  is  necessary  to  substitute  the  profile  given  by  equation  45  into  the  diffusion  differential 
equation  (one  for  each  constituent  of  the  atmosphere),  and  integrated  by  partial  fractions. 
This  was  done  by  Roberts  to  yield  the  density  for  the  altitude  band  between  100  and  125 
kilometers,  given  by  equation  52. 

P,(Z)  =  i>,(Z)  (52) 

iml 

It  is  computationally  expensive  to  calculate  the  density  at  100  kilometers  at 
different  exospheric  temperatures  by  using  equation  47.  Thus  values  for  the  density  at  100 
km  were  pre  computed  at  several  values  of  the  exospheric  temperature,  and  these  values 
could  then  be  extracted  instead  of  using  valuable  computer  time  and  memory.  Next  the 
atmospheric  constituent  mass  densities  are  calculated  by  the  following, 

P,(Z)=p0oO)^-m,  F"'‘ exp(M,tFj  (53) 

The  index  /  corresponds  to  the  values  1  through  6  of  the  various  atmospheric 
constituents  of  N2,  Ar,  He,  02,  O,  and  H,  respectively.  The  constants  found  in  equation 
53  are  the  characteristics  of  these  species  and  are  tabulated  in  Table  3  of  Lockheed  1992. 
Hydrogen  is  not  a  significant  constituent  below  125  kilometers,  so  it  is  not  included  in  the 


calculations.  (Lockheed,  1992,  p  B-89)  The  temperature  at  100  kilometers  is  calculated 
from  the  following, 


7tl00)  =  Tx  +Qdx  where  Cl  =  35~*  jrc„(l00)"  =  -0.94585589  (54) 


rp*Q 


and  F}  and  F4  are  calculated  by  equations  55  and  56 


^3  = 


z+K  1 

<?1 

r  Vi 

(  Z-',  1 

9l 

f  Z2-2XZ  +  X2+Y2  ) 

U*  +  iooJ 

o 

o 

1 

o 

o 

1 

[l002  -  200X  +  Jf24-f2J 

(55) 


75(z-ioo) 


+— tan' 


F *  {Z+ Ra){Ra  + 100)  Y 
where  the  parameters  q,  are  calculated  by 
1 


y(z-ioo) 


y2+{z-x){\oo-x) 


<12  = 


U{rx) 


(56) 


<?,  =  £/(',) 

l 

q.^l+r^Rl-X'-Y'  )q,  +  IV(r,  )q,  + 1 V(r2  )q, }  /  X’ 

<?6  =  ~<75  -  2(  X  +  Ra  )q„  -  (r2  +  F,  )<y,  -  (r,  +  Ra  )q2 
=  ~ 2^  ~q3~q2 

All  of  the  variables  in  these  equations  have  been  defined  earlier. 

For  the  region  above  1 25  kilometers,  it  is  still  a  valid  assumption  that  diffusive 
equilibrium  is  dominant,  but  the  temperature  given  by  equation  45  is  no  longer  valid. 
Jacchia  defined  the  temperature  profile  of  the  upper  altitude  regions  by  the  following 
empirical  asymptotic  function: 
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1 \Z)=TX  +5(r.-rr)tan-|o.95«||-;||^-SJl+45x  1 0 T(Z-  ] 2$“]|  (57) 

In  order  to  integrate  the  diffusion  differential  equations  in  closed  form,  Roberts  replaced 
equation  57  with  the  following  (Lockheed  1992,  p.  B-90) 


T\Z)  =  L-(Z-Tx)exp 


'  (Tx-TA(Z-\2S\ 
U-rJl  35  ) 


'  t  N 

\R*+ZJ 


(58) 


The  parameter  £  will  be  discussed  later. 

Integration  of  equation  58  yields  the  following  equation,  which  is  v«.Iid  for  all  of 
the  constituents  except  hydrogen. 


>,(Z)=p,  (125)^ 


T„-T 

Tm-Tx, 


(59) 


where 


r  _  M,gR< 

*  I 


R£T, 


( T~~TA(  35  ' 

Tx-T0  648 1.766, 


(60) 


At  this  point  several  corrections  are  made  for  the  particular  constituent  densities 
due  to  seasonal  changes.  The  value  of  the  helium  density  that  is  calculated  by  using 
equation  59,  must  be  corrected  due  to  the  seasonal  latitudinal  variation  as  given  by 
equation  44.  The  specific  form  is  presented  below 

L_,=P,(z)iola'-”’*  <6» 

where  /  corresponds  to  the  index  of  helium  presented  above.  Above  500  kilometers,  the 
concentration  becomes  significant,  therefore  it  must  be  accounted  for  by  the  following 


ft(Z)  =  p,  (500) 


”7t500)" 

l+or.+r. 

1 

1 

.  rtz) . 

71-7t500)J 

(62) 


where  the  hydrogen  density  at  500  kilometers  can  be  calculated  from 
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p(500)  =  ^W’,'M"4-!M’r~““,!'~l  (63) 

6  A 

The  temperature  at  500  kilometers  is  calculated  by  using  equation  58.  The 
constituents  are  summed  and  the  standard  density  above  125  kilometers  is  given  by  the 
following  (Lockheed,  1992,  p.  B-94): 

P,(Z)  =  5>,(Z)  (64) 

l»l 

So  far,  the  standard  atmospheric  density  at  any  given  altitude  has  been  calculated, 
but  the  densities  calculated  using  equations  47,  52,  or  64  must  now  be  corrected  for 
geomagnetic  activity,  the  semi-annual  variation  and  the  seasonal  latitudinal  variation  by 
equations  37,  38,  and  43  respectively.  The  effects  of  these  variations  can  be  summed 

logarithmically  to  obtain  the  following  relation 

(Alogp)^  =(Alogp)c  +(Alog  p)^  +(Alogp)tr  (65) 

The  final  corrected  density  at  any  altitude  can  then  be  determined  by 

p(Z)  =  pJ10(Alogp)~  (66) 

Due  to  the  introduction  of  equation  58  in  the  place  of  Jacchia's  equation  (57)  for 
temperature,  the  results  of  the  Roberts  portion  of  the  model  did  not  exactly  concur  with 
those  that  were  observed  by  Jacchia.  This  effect  was  only  in  the  upper  portion  of  the 
model,  whereas  the  lower  altitude  bands  agreed  exactly(as  should  be  the  case).  This  is 
where  the  (  parameter  comes  into  play  in  equation  58.  Values  of  the  C  parameter  were 
calculated  in  order  to  obtain  the  best  least  squares  fit  of  density  versus  altitude  (125  - 
2500km)  to  the  Jacchia  tabulated  data.  The  maximum  deviation  from  the  Jacchia  density 
values  is  less  that  6.7%.  (Lockheed,  1992,  p.  B-94)  The  derivation  and  values  of  the 
parameter  C  can  be  found  in  the  Lockheed  reference,  and  will  not  be  discussed  in  this 
paper. 
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a.  Mathematical  Basis 


The  following  section  contains  the  partial  derivatives  of  the  Jacchia- 
Roberts  model.  Starting  from  equation  66,  the  partial  derivative  of  density  with  respect  to 
local  position  can  be  found.  This  can  then  be  used  in  orbital  prediction  schemes  in  order 
to  better  model  the  satellite's  orbital  trajectory.  The  following  equation  is  a  simple 
restatement  of  equation  66. 

f(Z)  =  ps{Z) Ap_  (67) 


The  desired  partial  derivative  then  becomes 


9p 

dR 


=  Ps 


fapc) 

dR 


(68) 


The  variation  of  the  correction  factor  with  respect  to  the  satellite  position  is  derived  from 
equations  65  and  37  through  43. 


%hPc) 

dR 


Apc 


0.4342944819 


{*(/)  f(Z)S-+0  14sin(2ff0+ 1. 


(l-O.OO26{(Z-9O)}2}sin^jsin0)^+2(Z-9O)|sin^)cos0^ 

oR  oR 


(69) 


where 

/,(Z)  =  -O.002868/(Z)  +  2.33l(5.876xl0"7)z,33l^OO2868Z  (70) 

The  variation  of  altitude  with  respect  to  position,  dlj dR ,  is  the  same  as  that 
calculated  in  equation  24.  The  variation  of  geodetic  latitude  with  respect  to  position  is 
derived  by  differentiating  the  following  equation 


0  =  tan 


-i 


1 

*3 

(l-/)2 

1 

* 

K» 

+ 

& 

K> 

to  obtain 


(71) 


36 


dip  _  sin  20 
2 


- X , 


IT 


xr+x: 


x , 


*,2+xi 


1 


(72) 


The  variation  of  standard  density  with  respect  to  position  is  calculated  directly 
from  the  barometric  differential  equation  for  altitudes  below  100  kilometers 


dps 


— y  nA,Z-' 

DT 


RT 


dZ__}_dT_ 
dR  T  dR 


and  from  the  diffusion  differential  equation  above  100  kilometers 


(73) 


fa  =  il 

dR  T 


P'gpR 


R(Z+RJ 


dZ  ,  ,dT 

5K  +  p5  +  “3p3  dR 


where 


p'=±pM 

i=i 


(74) 


(75) 


The  partial  derivatives  of  the  temperature  with  respect  to  position  are  derived  by 
differentiating  equation  45  for  altitudes  less  than  125  kilometers 

(dTx  dT^ 


dT  T-  Tn 


dR 


dl :  dR 


+  1  A-i>CnZ"-'  1-pr 
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(76) 


or  equation  58  for  altitudes  above  125  kilometers 


dT  _  dl : 
dR  dR  + 


f  T-Tm 


1- 
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X  J 


dT 


Z-Zx 
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(  T- 71 

(71-?) 

«,  +  Zv 

1  <9Z 

'  35J 

'a/? 

(77) 


and  r  My,  the  derivatives  of  Tx  and  71  with  respect  to  position  are  derived  by 
ditic  .iitiating  equations  36  and  33  respectively, 

^  =  0. 05 1 8806 + (294. 3505)(0. 002 1 622)eM)002,6“  r-  (78) 

<971 

=  0. 3  Tc  L. :  Jn 1 1 2  0COS  1  -  cos3  0 
9R  c[  V  2  JdR 

12  •  30^  9r|  3  /  -» 2  -»  ■»  >  r  .  t  |  ._A. 

-2. 2 cos  rjsm  rjcos  — s-- (cos  q-sur*  0) cos*  -sin— -=4  (79) 

2  dR  2  2  2  ot?  J 
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3.  MSIS  86  (Mass  Spectrometer  and  Incoherent  Scatter  86) 

The  Lockheed  Densel  model  and  the  Jacchia  71  model,  both  developed  by 
Luigi  Jacchia,  use  drag  analysis  data  from  earth  orbiting  satellites  to  derive  the  densities  of 
the  various  altitude  bands  associated  within  the  models.  A.  E.  Hedin  chose  a  different 
approach  to  the  problem  of  modeling  the  atmosphere.  Hedin's  approach  was  to  use  actual 
density  data  retrieved  from  orbiting  instruments  aboard  several  satellites  and  combine  this 
data  with  measurements  taken  by  ground  based  incoherent  scatter  radar  stations.  Hedin, 
et.  al.,  began  the  original  MSIS  model  using  the  total  densities  inferred  from  Jacchia's 
models,  however,  it  was  found  that  the  measurements  of  the  temperature  found  by  the 
incoherent  scatter  method  differed  significantly.  (Hedin,  1972,  p.  2139)  These  differences 
were  noted  in  the  time  of  daily  maximum,  and  the  amplitude  of  the  daily  and  annual 
variations  of  the  density.  Hedin  used  the  data  that  was  obtained  from  the  mass 
spectrometers  on  board  several  orbiting  satellites  to  confirm  this.  The  measurements  from 
the  satellites  as  well  as  the  ground  based  incoherent  scatter  measurements  provide 
information  at  different  altitudes,  latitudes,  longitudes,  solar  activities ,  and  seasons. 

The  MSIS  86  model  is  the  culmination  of  years  of  research  and  development  as 
can  be  seen  in  figure  1.  The  formulation  of  the  MSIS  models  is  based  on  a  spherical 
harmonic  expansion  of  exospheric  temperature  and  effective  lower  boundary  densities 
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which  uses  local  solar  time  and  geographic  latitude  as  the  independent  variables.  This 
expansion  is  a  special  case  of  a  more  general  expansion  based  on  longitude  and  latitude, 
where  the  coefficients  of  the  expansion  are  represented  by  a  Fourier  series  in  universal 
time.  (Hedin,  1979,  p.  2) 

This  leads  to  the  following  expansion 

AG  =  XX5>.  [a£  co^ma)'t+nX)  +  b”  sin(wa>7+//A)]  (80) 

npmO  ImQ  #»*-/ 

where 

Pto  =  Legendres  associated  functions  in  geographic  latitude 
t  -  universal  time  in  seconds 
m,  =  2ff/86400j-1 

X  =  geographic  longitude  in  radians 

This  particular  expansion  was  used  for  the  initial  MSIS  models.  With  the 
launching  of  several  more  satellites  equipped  with  mass  spectrometers,  and  the  addition  of 
several  more  incoherent  scatter  radars,  the  model  evolved  into  the  MSIS83,  and  eventually 
the  MSIS  86. 

The  main  emphasis  now  is  the  formulation  of  the  MSIS  86  model.  Hedin 
explains  that  the  models  uses  a  Bates  temperature  profile  as  a  function  of  geopotential 
height  for  the  upper  thermosphere,  and  an  inverse  polynomial  in  geopotential  height  for 
the  lower  thermosphere.  (Hedin,  1987,  p.  4649)  The  exospheric  temperature  and  other 
key  quantities  are  expressed  as  functions  of  geographical  and  soiar/magnetic  parameters. 
The  temperature  profiles  allow  for  the  exact  integration  of  the  hydrostatic  equation  for  a 
constant  mass  to  determine  the  density  profile  based  on  a  density  specified  at  120 
kilometers  as  a  function  of  geographical  latitude  and  solar/magnetic  parameters.  The 
MSIS  83  model  used  the  expansion  formula  (Equation  80)  to  model  the  variations  due  to 
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local  time,  latitude,  longitude,  universal  time,  FI 0.7,  and  Ap.  The  MSIS  86  model 
enhances  the  MSIS  83  model  by  adding  terms  to  express  hemispherical  and  seasonal 
differences  in  the  polar  regions  and  local  time  variations  in  the  magnetic  activity  effect. 

The  following  equations  are  the  complete  formulation  of  the  MSIS  86 
atmospheric  model.  The  first  step  in  the  formulation  is  to  obtain  an  expression  for  the 


temperature  profile. 

7t  Z)  =  £-(£- 7;)  exp^22''1  where  Z>Zfl  (81) 

7{Z)  =  l/(l/r0  +  Tbx2  +  Tcx 4  +  Tdx6  )  where  Z  <  Z„  (82) 

the  matching  temperature  and  temperature  gradient  at  Za  equals 

T„  =  T[Zj  =  T.  -(r.  -  ^exp^12-*’1  (83) 

t’  =  r(Zj=(T„  -  Tj^(Rr  +  Zl)/(Rr+Z,)]‘  (84) 

Ta  =  0.66666£(Zo,Zj  C/^2 -31111  l(l/7;  -  1/fJ  +  7. 1 11 1 1(1/7;,  - l/T„)  (85) 

Tc  =  4(Z„ZjT:/(2Tj‘)-(\/Z-\/n)-2TD  (86) 

T,=(l/T'-n)-rc-T0  (87) 

x  =  -[fe.Z)-  KZ„Zj]/|(Z0,Zj  (88) 

T„=T,  +  Tr(T,-T,)  (89) 

KZ,Z,)=(Z-Z,)(«,+Z,)/(7?,+z)  (90) 

KZ.Zj  =  (Z-Zj(«,  +  zJ/(/!,+z)  (91) 

/(t„-t,)  t;  =  77[i+G(jT.)} 

j  =  t'[U-G(L)}  7"0  =  5[l+G(L)] 

7.  =  r[l+G(£)]  Z,  =  Z,[l  +  G(l)] 

r,  =  fs[i+G(£)] 

where  (all  temperatures  in  Kelvin) 


T = ambient  temperature  T0  =  average  mesopause  temperature 
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Ta  =  temperature  at  Za  T,  =  average  temperature  gradient  at  Z, 

i  ___ 

Ta  =temperature  gradient  at  Za  in  K/km  TR  =  average  mesopause  shape  factor 
Tn  =  temperature  at  (Z0  +  Z, )  /  2  Z  =  altitude  in  km 

Tm  =  average  exospheric  temperature  Z,  =  120  km 
Tt  =  average  temperature  at  Z, 

Za  =  altitude  of  temperature  profile  junction;  1 17.2  km 

Z0  =  average  mesopause  height  Rp  =  6356.77  km 


Hedin  now  explains  that  the  density  profile  is  a  combination  of  diffusive  and 
mixing  profiles  multiplied  by  one  or  more  factors  C, . . .  C„  to  account  for  chemistry  or  the 
dynamics  flow  effects.  (Hedin,  1987,  p.  4638) 

;KZ,A/)  =  [»,(Z.A/)''+«m(Z,A/)'l]','*[C,(Z)...C.(z)]  (92) 

A  =  (93) 

where 


n  =  ambient  density 
nd  =  diffusive  profile 
nm  =  mixing  profile 


Mi  =  28 

M0  =  28.95 

M  =  molecular  weight  of  gas  species 


The  Following  equations  represent  the  diffusive  profile  and  the  mixing  profile. 
The  diffusive  profile  is  given  by  the  following 

na(Z,M)=n,D(Z,M)[l\Z,)ll\Z)\'°  (94) 

D(Z,M)  =  D,(Z,M)  ZSZ.  (95) 

d(z,m)= d,(z. , (z<z)  {96) 

D,(ZM)  =  \T(Z,)ll{z)f  expl*w•,<z•2,,|  (97) 
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>1 

II 

ga=g,/(l  +  Za/Rp)' 

Yi=Mga${Z0,Za)/Rg 

g,  =g,/(i+Zt/RPY 

n,  =  fi,  exp;G,z^ 

(98) 

where 

77,  =  the  average  density  at  Z, 

gt  =9.80665x10 -'km/s2 

Rg  =  8. 3 14  x  10"3  g  km^/mols  deg 

and  a  is  the  thermal  diffusion  coefficient  of  -0.4  for  He  and  zero  for  the  rest  of  the 
constituents. 

The  mixing  profile  is  given  by  the  following 
».(Z,M)  =  »,[z>(Z,>M)/o(z,J 70)|r(Z,)/7tZj]"fl(Z,iK,)[7(Z,)/7tZ)]  (99) 


Hedin  (1987),  goes  into  a  procedure  to  simulate  the  chemistry  and  dynamic  flow 
effects  of  the  turbopause  which  provides  a  specified  mixing  ratio  with  respect  to  nitrogen. 

C,(Z)  =  exp{/?,^l  +  exp{z~z,)/H'  jj  (100) 

fi,  =log[Q«.(Z„28)/«.(Z„W)]  (101) 

where  £1  is  the  mixing  ratio  relative  to  N2,  Z  j  is  the  altitude  at  which  log  10  Cj  is  R/2, 
and  Hi  is  the  scale  height  for  this  correction.  There  is  a  peak  in  the  mixing  ratio  in  the 
lower  thermosphere  due  to  O,  N,  and  H.  This  peak  is  modeled  by  the  following 

C:(z)  =  exp{«,/[l  +  exp(Z-Z!)///!]}  (102) 

where  R2  is  the  density  correction  parameter,  Z2  is  the  altitude  where  log  10  density 
correction  is  R2/2,  and  H2  is  the  scale  height  of  the  correction.  As  mentioned  above,  the 
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MSIS  86  atmospheric  model  attempts  to  model  all  the  known  causes  of  the  variations  of 
density.  This  is  done  by  the  expansion  function  (Equation  80).  The  following  is  the 
expansion  function  for  the  MSIS  86  model  quantities. 

Time  independent 

G  =  Q\qP\q  +  fl;o4  ■*"<I40^40 

Solar  activity 

+F£6F  +  F£(t>Fy  +f£AF+/£(bF)1  +/;'/>„  AF 

Symmetrical  annual 

+ci,  cosfl,  (/,-/£) 

Symmetrical  semiannual 

+(4  +44) C0S2O,  (td  -  4* ) 

Asymmetrical  annual  (seasonal) 

+(cio4  +  c3o4)^*i  cos2£2</(frf  ~4) 

Asymmetrical  semiannual 

+cio4  cos2 Qd{ld  ~*w ) 

Diurnal 

+[ai]4  +  fl3l4  +  fl5l4  +  C2l4  COS^j(0  ^10  C0S^ 

+[44  +  44  +  44  +  4l4  C0S^d(^d  — ^to)]^2  Sm  tOT 

Semidiurnal 

+[a22  4  +a424  +(44  +44)cos  £ij(td  —  4o)]4  cos2cor 

+[^22 4  +44  +(^32 4  +  £44)C0S^</(*«/  ~  ^1*0  )]4  sin2tur 
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Terdiumal 


+[a33^33  +(c'*Ai  +^iP6i)cosQ.J{tJ  ~^)}fz  cos 3 cur 
"*■[^33^33  +  (^43^43  +  ^63^63  )  COS^ij  (tj  ~  f|o  )]^*:  S*n  3CUT 
Magnetic  activity 

+{*i  +k^Px+k’„P„  +(*J3>,  +  k£P„  +</>10)cos£iJ(/J  -C‘) 


Longitudinal 


h[<4/>;,  +<3>,  +<.“/>„  +<.»/>,  +a»f!1+a«p>l  +«/;,  +<3,„)co»QJ(/<-/S)] 
x(l+F,10AF)cosA 


+[«:",a, +»x  +*x  +*5°,-p3,  +ta  +»£/’„  )cosii„(<,  -/,■;)] 

x(l  +  F2“0AF)sinA 
UT 

+te+<P„+o^„)(i+^'Af)(i+rl;^J(i+r,;,cosnJ(/J-C)) 

XCOSO)'((  -  C  )  +(<4  P„  +«!;  4  )  COi[<l»'('  -  ^  )  +  2A] 

UT/longitude/magnetic  activity 

+(k"P„  +*;Pn  +C^JO+^XK<cos(A-A*°) 

+*,‘X  cosQ,,(i,,  -  rfd  )AA  cos(A  -  A“ ) 

+(*,;' C  +  *£/>„  +  *i'  />„  ]  A4  cos<B-(<  -  /“ )  (103) 


f|=l+^AF+/,;,AF+/(;!(AF)i  (104) 
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where 


F2=l  +  Fl‘AF+f£AF+/‘2(AF)2 


(105) 


AF  =  FW1-Fl01  AF  =  Fl0J-l5Q  F10.7  is 

the  solar  flux  measurement  for  the  previous  day,  F]01  is  the  81  day  average  of  the  solar  flux 
measurement  centered  on  the  day  in  question. 

Pnm  are  the  non  normalized  Legendre-associated  functions  equal  to  the  following 

[(l-r2)*2  l2Hn\^dn+mldxn+m)(x2  - 1)" 

with  x  equal  to  the  cosine  of  the  geographic  latitude, 

Qd  =2tt/365  (o  —  2jt/24  co'  =  2^/86400,  x  is  local  time  in  hours,  t  is 
universal  time  in  seconds,  t<j  is  the  day  count  in  year,  and  X  is  the  geographic  longitude 
with  eastward  being  positive.  Hedin  explains  that  there  are  two  choices  for  magnetic 
activity, 

AA  =  ( Ap  - 4) + - 1){ Ap  -  4 + [exp( -k‘n ( Ap  - 4))  - 1] / }  ( 1 06) 

where  Ap  is  the  daily  magnetic  index,  or  there  is  another  method  illustrated  that 
establishes  the  magnetic  activity  that  uses  the  average  values  over  an  extended  period  of 
time.  (Hedin,  1987,  p.  4661) 
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III.  ATMOSPHERIC  MODEL  EVALUATION 


The  evaluation  of  atmospheric  drag  models  is  a  relatively  difficult  procedure  owing 
to  the  fact  that  there  is  no  real  comparison  tool  to  be  used.  Until  enough  satellites  can  be 
placed  into  various  orbits  and  each  of  these  satellites  carries  density  measuring  equipment, 
the  density  of  a  particular  orbit  cannot  truly  be  obtained.  As  mentioned  previously,  the 
modeling  of  the  earth's  atmosphere  is  the  most  inaccurate  of  all  the  perturbations 
encountered  in  orbital  motion  analysis.  The  intent  of  this  project  was  to  obtain  knoton 
vectors  of  a  particular  satellite  program,  and  run  a  comparison  of  an  orbital  prediction 
scheme  using  each  of  three  atmospheric  models  for  two  different  altitudes.  In  order  to 
make  the  comparison  valid,  consistency  of  the  prediction  scheme  had  to  be  verified. 

The  vectors  being  used  in  the  analysis  were  obtained  from  the  Air  Force  Satellite 
Tracking  Center  on  two  satellites,  one  at  650  kilometers  and  one  at  approximately  800 
kilometers,  both  of  which  have  a  high  inclination.  The  vectors,  position  and  velocity 
represented  in  True  of  Date  Cartesian  coordinates,  are  obtained  by  triangulating  the 
satellite  position  and  predicting  out  to  20:00:00.00  local  time.  It  would  have  been  better 
for  the  evaluation,  if  the  actual  observed  position  of  the  satellites  could  have  been 
obtained,  but  due  to  time  constraints  and  logistical  problems  this  could  not  be  achieved. 
The  accuracy  of  the  position  and  velocity  vectors  is  advertised  by  the  STC  to  be  within 
2000  ft  in  the  X,  Y,  and  Z  directions.  It  was  decided  at  the  time,  that  this  accuracy  was 
acceptable  for  the  evaluation. 

In  order  to  provide  consistency,  the  prediction  scheme  developed  had  to  be  the  same 
as  that  used  by  the  STC.  This  proved  to  be  quite  a  troublesome  problem.  The  STC  uses  a 
prediction  scheme  which  models  most  if  not  all  of  the  known  perturbations.  Depending  on 
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what  is  requested,  certain  perturbations  can  be  temporarily  removed  from  the  program. 
The  vectors  obtained  from  the  STC  contained  only  those  perturbations  which  were 
inherent  to  the  prediction  scheme  developed,  with  one  exception.  As  mentioned 
previously,  the  STC  uses  the  WGS-84  geopotential  model.  A  great  deal  of  effort  was 
placed  on  placing  the  same  geopotential  model  in  the  developed  prediction  scheme,  but 
several  problems  soon  became  apparent.  The  conversion  of  True  of  Date  coordinates  into 
the  WGS-84  frame  of  reference  is  not  a  trivial  process.  An  attempt  was  made  to 
,  incorporate  the  full  41  by  41  geopotential  matrix  in  the  scheme,  but  error  build  up  and 
computational  time  increased  greatly.  It  was  decided  to  incorporate  a  truncated  version  of 
the  WGS-84  geopotential  model.  A  six  by  six  matrix  representing  the  first  six  zonals  of 
the  model  was  eventually  used  to  prevent  error  build  up.  Future  work  on  this  project 
would  have  to  include  the  conversion  of  the  TOD  coordinates  to  the  WGS  frame  of 
reference  and  incorporation  of  the  full  41  by  41  matrix. 

Another  way  that  consistency  was  achieved  was  to  verify  and  use  astrodynamical 
constants  consistent  with  AIAA  standards.  Depending  on  which  reference  one  uses  to 
obtain  values  such  as  the  earth's  gravitational  parameter,  earth's  equatorial  radius,  earth 
flattening  constant  and  others,  these  values  would  all  vary  by  some  small  amount.  It  was 
found  that  depending  on  which  constant  was  used,  the  error  would  be  different.  It  was 
decided  to  use  the  AIAA  published  values  whenever  possible.  These  values  were  placed 
in  a  subroutine  of  the  prediction  scheme  to  be  used  as  necessary. 

Once  it  had  been  determined  that  the  developed  prediction  scheme  was  an  equal 
representation  of  the  STC's  prediction  scheme,  or  as  close  as  it  could  be  gotten  in  the  time 
allowed,  the  model  evaluation  was  begun.  The  vectors  obtained  from  the  STC  spanned 
the  time  frame  from  21  August  1993  to  22  September  1993  in  24  hour  forecasts.  That  is, 
the  position  and  velocity  were  reported  for  20:00:00.00  of  each  day.  These  vectors  were 
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then  programmed  into  60  satellite  files,  30  for  the  650  kilometer  satellite  and  30  for  the 
800  kilometer  satellite.  Each  of  the  satellite  files  contains  not  only  the  position  and 
velocity  vectors,  but  also  the  environmental  conditions  of  FI  0.7  and  Ap  for  that  day,  plus 
the  8 1  day  average  FI 0.7  as  reported  by  the  National  Geophysical  Data  Center.  The 
vectors  were  propagated  over  a  24  hour  period  and  compared  to  the  next  day's  vector. 
Since  the  baseline  vectors  themselves  contained  error,  it  was  decided  to  compare  the 
results  by  using  the  specific  mechanical  energy,  a  conserved  quantity,  of  each  of  the 
satellites  over  the  30  day  period.  By  converting  the  position  and  velocity  of  the  STC 
vectors  and  those  obtained  from  the  developed  prediction  scheme  into  specific  mechanical 
energy  using  equation  107,  a  comparison  could  be  made. 


(107) 


A  relative  comparison  of  the  position,  velocity,  right  ascension,  declination,  and  radial 
distance  was  then  conducted  in  order  to  determine  which  model  provided  the  most 
accurate  results,  the  object  being  to  meet  or  beat  the  2000  foot  accuracy  advertised  by  the 
STC  prediction  scheme. 

As  mentioned,  the  evaluation  procedure  consisted  of  propagating  the  baseline  vector 
for  a  24  hour  integration  period  for  each  of  the  satellites,  then  switching  atmospheres  and 
running  the  routine  again.  It  must  be  noted  that  the  Jacchia  60  and  Jacchia  71  models 
were  run  with  relatively  little  difficulty.  The  MSIS  86  atmospheric  model  proved  to  be 
quite  the  opposite.  Due  to  an  unexpected  hardware  failure  late  in  the  project,  and  the  fact 
that  the  atmospheric  model  itself  is  quite  difficult  to  integrate  into  the  prediction  scheme, 
the  evaluation  of  this  model  could  not  be  completed  in  time.  Work  is  currently  under  way 
to  correct  this  situation. 
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A.  JACCHIA  60  MODEL  EVALUATION 

The  first  evaluation  was  conducted  on  the  Jacchia  60  earth  atmospheric  model.  The 
first  procedure  was  to  determine  the  values  for  the  solar  flux  and  geomagnetic  activity 
encountered  during  the  time  period  between  August  21  to  September  22,  1993.  Figure  13 
shows  the  reported  values  of  FI  0.7  and  Ap  for  this  time  period.  As  compared  to  Figure  6, 
it  can  be  seen  that  the  solar  flux  values  are  on  the  low  side  of  the  solar  cycle,  and  the 
geomagnetic  activity  is  relatively  low.  This  allows  for  an  easier  comparison,  because  a 
perturbed  atmosphere  would  have  made  some  of  the  data  points  obtained  during  an  active 
period  out  of  range  of  those  obtained  during  a  quiet  period. 

The  first  evaluation  was  conducted  on  the  satellite  at  650  kilometers.  Figure  14 
illustrates  the  specific  mechanical  energy  of  the  STC  reported  vectors  versus  the  predicted 
vectors.  It  can  be  seen  that  the  predicted  vectors  are  extremely  close  to  that  of  the  STC 
reported  values.  The  two  error  bars  located  on  the  graph  illustrate  the  2000  ft  error  range 
advertised  by  the  STC  prediction  scheme.  In  one  case  it  was  found  that  the  STC  obtained 
vector  on  Julian  date  242  was  actually  reported  incorrectly.  For  the  most  part,  however, 
the  predicted  values  showed  only  small  variations  with  the  baseline  vectors.  In-track 
errors  fell  within  the  range  of  1700  to  9000  feet  in  radial  distance,  and  0.01  to  0. 15 
degrees  in  declination.  Cross-track  errors  also  fell  within  the  range  of  0.01  to  0. 15 
degrees.  Satellite  velocities  had  some  small  error  in  the  range  of  0.2  to  ten  feet  per 
second. 

The  same  procedure  was  then  conducted  on  the  800  kilometer  satellite  vectors. 

Once  again  the  prediction  scheme  provided  vectors  very  similar  to  those  obtained  from  the 
STC.  Figure  15  illustrates  the  results.  It  must  be  noted  at  this  point,  that  the  accuracy  of 
the  prediction  scheme  suffered  at  the  higher  altitude.  This  is  due  to  the  atmospheric 
model.  As  mentioned  previously,  accuracy  of  the  atmospheric  model  deteriorates  the 
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higher  the  altitude.  In-track  errors  were  found  to  range  from  2500  to  9000  feet  in  radial 
distance,  and  0.01  to  0.2  degrees  in  declination.  Cross-track  errors  were  in  the  range  of 
0.01  to  0.8  degrees.  Velocity  was  in  the  range  of  one  to  eight  feet  per  second. 

B.  JACCHIA  71  MODEL  EVALUATION 

Once  again  the  procedure  was  run,  this  time  with  the  Jacchia  7 1  model  inserted  into  the 
prediction  scheme.  Figure  16  illustrates  the  results  of  the  prediction  runs  on  the  650  kilometer 
satellite.  The  prediction  scheme  provided  very  little  variation  with  the  baseline  vectors.  The  in¬ 
track  errors  were  250  to  8000  feet  in  radial  distance  and  0.001  to  0. 14  degrees  in  declination. 
Cross-track  errors  were  also  0.001  to  0. 14  degrees  in  right  ascension.  Errors  in  velocity  predicted 
were  0. 1  to  ten  feet  per  second.  The  results  of  this  atmosphere  at  650  kilometers  were  much  more 
accurate  than  those  of  the  Jacchia  60  model.  In  order  to  confirm  this,  the  prediction  scheme  was 
run  again  with  the  vectors  of  the  850  kilometer  satellite.  In-track  errors  were  in  the  range  of  200  to 
8000  feet  in  radial  range  and  0.001  to  0.10  degrees  in  declination.  Cross-track  errors  were  also 
comparable  with  the  error  range  being  0.001  to  0. 10  degrees  of  right  ascension.  This  confirmed 
the  fact  that  the  Jacchia  7 1  model  provided  better  accuracy  at  both  altitudes.  This  is  due  to  the  fact 
that  the  Jacchia  71  model  provides  a  more  accurate  modeling  of  the  polar  regions,  whereas  the 
Jacchia  60  does  not  model  the  polar  regions  as  well.  The  modeling  of  the  polar  regions  is 
extremely  important  for  the  satellite  program  being  studied  due  to  its  high  inclination. 
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IV.  SUMMARY 


This  thesis  has  attempted  to  provide  the  underlying  principles  used  in  orbital 
prediction  and  a  comparison  of  atmospheric  models.  The  prediction  scheme  that  has  been 
created  was  used  to  compare  the  Jacchia  60  and  Jacchia  71  earth  atmospheric  models.  It 
was  found  that  the  Jacchia  71  model  provided  much  more  accurate  results  than  the  Jacchia 
60.  This  would  stand  to  reason  due  to  the  fact  that  the  J71  model  had  a  bigger  data  base 
with  which  to  work  from.  The  J71  model  provides  density  information  that  is  dependent 
both  on  solar  activity  as  well  as  geomagnetic  activity,  whereas  the  J60  only  relies  on  solar 
activity.  The  J71  also  models  the  polar  regions,  which  is  essential  to  the  satellite  program 
being  investigated.  The  MSIS  86  earth  atmospheric  model  could  not  be  compared  due  to 
time  constraints,  but  is  currently  being  integrated  into  the  prediction  scheme. 

Several  improvements  can  be  made  to  the  prediction  scheme  in  order  to  improve 
accuracy.  As  discussed  earlier,  the  geopotential  model  needs  to  be  incorporated.  The 
True  of  Date  coordinates  need  to  be  transferred  into  the  WGS-84  frame  of  reference,  and 
the  complete  41  by  41  geopotential  matrix  incorporated.  Doing  this  would  increase 
computational  time,  but  accuracy  would  greatly  improve.  Another  improvement  in  the 
prediction  scheme  would  be  to  vary  the  B-factor  as  the  satellite  cross-sectional  area 
changes  through  its  orbital  pass.  Currently  the  B-factor  is  held  at  a  constant  value.  By 
changing  the  B-factor  as  a  function  of  latitude  or  declination,  error  can  be  reduced. 

Follow  on  work  on  this  project  would  include  inserting  other  atmospheric  models 
for  comparison,  such  as  the  MSIS  86  and  MSIS  90.  Before  doing  this,  however,  it  would 
be  better  if  the  satellite  vectors  being  used  were  actual  observations  and  not  predictions 
themselves.  This  would  provide  a  better  realization  of  the  accuracy  of  the  program,  and 
hence  the  atmospheric  models.  Another  improvement  in  the  program's  accuracy  would  be 
to  convert  from  True  of  Late  Cartesian  coordinates  to  normalized  spherical  coordinates 
and  then  integrate.  Doing  this  would  greatly  decrease  the  error  build  up.  As  the 
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prediction  scheme  stands  now,  the  error  build-up  causes  accuracy  problems  if  the  vectors 
are  integrated  for  more  than  a  twenty  four  hour  period. 

Orbital  prediction  is  an  essential  tool  in  today's  space  age.  in  the  case  of  the  satellite 
program  of  interest,  if  the  accuracy  of  the  orbital  prediction  can  be  improved,  then  the 
mission  analysis  can  be  improved.  The  interesting  point  of  this  project,  was  that  the  B- 
factor  was  no  longer  being  used  as  the  error  catch-all.  The  B-factor  for  any  given  satellite 
can  be  correctly  used,  and  the  position  and  velocity  as  well  as  any  other  of  the  satellite 
ephemirides  can  be  calculated.  By  incorporating  more  subroutines  that  model  other 
perturbing  forces,  this  prediction  scheme  can  be  used  for  other  satellite  programs. 
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Figure  1.  Atmosphere  Breakdown 
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Figure  2.  Temperature  vs.  Altitude 
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Figure  3.  Temperature  Profile  of  Earth's  Atmosphere 
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Figure  4.  Height  Variation  of  Number  Density 
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Figure  6.  Solar  Flux  History 
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Figure  8.  Solar  Radiation  Spectrum 
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Figure  9.  Geomagnetic  Storm  Recording 
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Figure  11.  Comparison  of  Perturbation  Magnitudes 
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Figure  12.  Atmospheric  Model  History 
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Daily  Measured  Values  of  FI  0.7cm  and  Ap 


Figure  13.  Observed  values  of  FI  0.7  and  Ap  for  21  August  to  22  September  1993 
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Figure  14.  Jacchia  60  Evaluation  at  650  km 
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Figure  17.  Jacchia  71  Evaluation  at  800  km 
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TABLE  2:  TYPICAL  BALLISTIC  COEFFICIENTS  FOR  LEO  SATELLITES 


Satellite 

I 

Shape 

Max. 

cross- 

sectio¬ 

nal 

arts 

(ma) 

Min. 

cross- 

sectio 

-net 

area 

(mJ> 

Drag 

coettl- 

ciant 

lor 

max. 

cross. 

arts 

Drag 

coeffi¬ 

cient 

lot 

min. 

cross. 

area 

Max. 

ballistic 

coeffi¬ 

cient 

(kgim1) 

Min. 

ballistic 

coeffi¬ 

cient 

(kglm9) 

Type  of 
Mission 

Oscar- 1 

5 

box 

0  075 

00584 

4 

2 

428 

167 

comm 

Intercosmos 

-16 

550 

cylin 

2.7 

3.16 

2.67 

2  1 

82  9 

76  3 

scientific 

Viking 

277 

octag 

225 

0833 

4 

26 

128 

308 

scientific 

Explorer- 11 

37 

octag 

0.18 

007 

283 

26 

203 

72  6 

astronomy 

Explorer- 17 

1882 

sphere 

0621 

0621 

2 

2 

152 

152 

scientific 

Space  Tele¬ 
scope 

11000 

cylmd.  w/ 
arrays 

112 

143 

333 

■ 

192 

295 

astronomy 

OSO-7 

834 

9-sided 

105 

05 

367 

29 

437 

165 

solar  phys. 

OSO-8 

1063 

cylmd  w / 
arrays 

599 

1  81 

376 

■ 

147 

472 

solar  phys 

Pegasus-3 

10500 

cylind  w/ 
arrays 

264 

14  5 

33 

■ 

181 

12  1 

scientific 

Landsat-1 

891 

cylind.  w/ 
arrays 

10  4 

1  81 

34 

■ 

123 

252 

remote 

sensing 

EBS-1 

2160 

box  w  / 
arrays 

45  1 

m 
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■ 

135 

120 

remote 

sensing 

LOEF-1 

9695 

39 

14  3 

267 

n 

169 

93  1 

exper 

platform 

HEAO-2 

3150 

fiihiB'B 

■ 

139 

4  52 

2  83 

Bi 

174 

80  1 

astronomy 

Vanguard-2 

939 

sphere 

0.2 

2 

2 

235 

235 

scientific 

SkyLab 

76136 

cylind  w/ 
arrays 

462 

464 

35 

■ 

410 

47  1 

scientific 

Echo-1 

75  3 
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731 

731 

2 

2 

0  515 

0.515 

comm 

Extrema 

437 

0  515 
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